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A  new  approach  for  assessing  the  survivability  of  air¬ 
craft  components  in  nnclear  blast  and  thermal  environments 
is  presented  in  this  dissertation.  A  nonparase trie  tech¬ 
nique  is  developed  for  use  in  calculating  the  reliability 
interference  integral.  This  approach  eliminates  the  need 
for  density  function  identification  and  parameter  estima¬ 
tion.  Furthermore,  the  method  can  be  used  without  resorting 
to  large-sample  random  Monte  Carlo  simulation  or  propagation 
of  moments.  In  addition  to  this,  the  derived  cumulative 
distribution  function  using  snch  a  technique  exactly  inter¬ 
polates  the  true  distribution  function  at  selected  points. 
The  method  is  applied  to  the  problem  of  aircraft  survivabil¬ 
ity  in  nuclear  blast  environments  using  failure  (strength) 
distributions  found  in  the  literature.  It  is  also  applied 
to  the  case  of  aircraft  survivability  in  nuclear  thermal 
environments  where  direct  failure  data  is  not  available. 
Inputs  to  the  engineering  models  involved  are  treated  sta¬ 
tistically,  and  the  method  is  used  to  rigorously  determine 
the  statistical  nature  of  the  output  variables. 
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The  objective  of  this  dissertstion  is  to  find  a  logical 
way  of  inferring  the  probability  of  failure  of  an  aircraft 
component  when  exposed  to  the  effects  of  nuclear  weapons.  A 
reliability  theory  approach  is  taken  that  is  distribution 
free.  This  approach,  which  exploits  some  recent  develop¬ 
ments  in  nonparametr ic  statistics,  eliminates  the  need  for 
density  function  identification  and  parameter  estimation. 
Furthermore,  it  provides  a  new  way  of  finding  the  distribu¬ 
tion  of  a  function  of  random  variables  without  using  large- 
sample  Monte  Carlo  simulation.  The  method  is  applied  to  the 
problem  of  aircraft  survivability  in  nuclear  blast  environ¬ 
ments  using  test  data  found  in  the  literature.  It  is  also 
applied  to  the  case  of  aircraft  survivability  in  nuclear 
thermal  environments  where  direct  failure  data  is  not  avail¬ 
able.  Inputs  to  the  engineering  models  involved  are  treated 
statistically,  and  the  method  is  used  to  rigorously  deter¬ 
mine  the  statistical  nature  of  the  output  variables. 

A  review  of  nuclear  s ur v i v a b i 1  i  ty / vulner ab il i  ty  (S/V) 
methods  is  provided  in  Chapter  II.  The  deterministic 
approach  for  aircraft  is  reviewed.  Little  work  on  probabil¬ 
istic  methods  for  aircraft  has  been  done  up  to  this  time. 
Attempts  to  logically  develop  a  probabilistic  model  from 
stated  aircraft  hardness  specifications  are  noted  as  are 
other  approaches.  Ground  system  survivability  methods  are 


examined.  Probabilistic  methods  have  been  used  for  some 
time  in  this  area.  The  earliest  approaches  involve  direct 
modeling  of  system  failure  probability  as  a  function  of 
range  from  the  weapon.  Modern  methods  include  Monte  Carlo 
simulation.  This  has  been  used  in  Minuteman  survivability 
studies  and  more  recently  to  investigate  electromagnetic 
pulse  (EMP)  survivability.  Other  studies  include  surviva¬ 
bility  assessment  of  horizontal  and  vertical  missile  shel¬ 
ters  and  interruption  of  low  frequency  communications  by 
nuc 1  ear  effects. 

A  review  of  classical  s t r e s s / s t r e ng t h  interference 
theory  is  provided  in  Chapter  III.  Strength  and  stress 
distributions  are  discussed.  and  the  reliability  interfer¬ 
ence  integral  presented.  Engineering  determinism.  or 
"cookie-cutter'’  damage  distributions,  are  shown  to  be  a 
special  case  of  the  general  interference  theory  approach. 
Stress/strength  interference  theory  is  thus  a  natural  choice 
as  a  theoretical  approach.  Serious  problems  have  existed  in 
applying  interference  theory  to  failure  assessment  of  large 
engineering  systems.  These  include  development  of  the  sys¬ 
tem  reliability  model  from  component  structures,  the  numeri¬ 
cal  difficulty  of  finding  the  distribution  of  a  function  of 
one  or  several  random  variables,  and  the  serious  limitation 
presented  by  the  lack  of  data.  This  last  limitation  in  par¬ 
ticular  makes  density  function  identification  and  parameter 
estimation  difficult.  Attempts  to  solve  the  above  diffi- 


culties  include  fault  tree  analysis,  propagation  of  sonant 
methods,  direct  and  indirect  Konte  Carlo  simulation,  vari¬ 
able  transformation  techniques,  Bayesian  analysis,  and  sur¬ 
vey  of  expert  opinion. 

A  new  solution  to  the  above  difficulties  using  an 
extension  of  recently  developed  no np a r as e t r i c  estimation 
techniques  is  presented  in  Chapter  IV.  While  the  systes 
reliability  modeling  problem  remains,  the  new  numerical 
technique  provides  a  tool  for  finding  the  distribution  of  a 
function  of  multiple  random  variables.  An  example  from  the 
engineering  literature  is  presented  as  a  benchmark  problem. 
The  advantages  of  the  nonparametric  method  are  noted.  These 
advantages  include  (a)  elimination  of  the  requirement  for 
density  function  identification  and  consequent  parameter 
estimation,  (b)  freedom  from  a  priori  distributional  assump¬ 
tions,  (c)  elimination  of  the  requirement  for  large  Monte 
Carlo  samples,  and  (dj  protection  against  drawing  unwar¬ 
ranted  inferences  from  a  small  data  set. 

The  analysis  of  aircraft  survivability  to  nuclear 
induced  blast  environments  is  presented  in  Chapter  V.  In 
this  case,  actual  failure  distributions  for  piece-parts  are 
taken  from  an  analysis  of  some  20  years  of  static  test  data 
conducted  at  Wright-Patterson  AFB,  OB.  A  general  discussion 


tail  assemblies.  This  is  done  by  nsing  the  statistical 
▼ariation  of  overpressure  versus  range  to  determine  the 
applied  stress  distributions.  The  actual  failure  distribu¬ 
tions  determined  from  the  test  data  then  allow  the  component 
failure  probabilities  to  be  determined  as  a  function  of 
range.  These  provide  the  required  information  for  finding 
the  system  failure  probability  as  a  function  of  range  from 
the  weapon. 

The  analysis  of  aircraft  survivability  to  nuclear 
thermal  environments  is  presented  in  Chapter  VI.  The  diffi¬ 
culty  of  sparse  data  is  noted.  In  particular,  no 
probabilistic  description  of  the  thermal  environment  could 
be  found  from  the  literature.  Also,  no  failure  data  for  the 
failure  mode  of  interest  could  be  found  in  the  literature. 
As  an  alternative,  the  statistical  distribution  of  radiated 
power  is  developed  by  regressing  the  output  of  a  complex 
radiation  hydrodynamics  code  on  that  of  a  more  approximate 
model.  The  study  is  restricted  to  the  radiated  power  from  a 
1  Kiloton  (KT)  burst  in  a  sea  level  atmosphere.  The  statis¬ 
tical  environmental  input  is  then  used  to  determine  the 
distribution  of  peak  shin  temperature  in  a  thin  shin  assem¬ 
bly.  Two  thermal  failure  models  are  compared,  and  failure 
probabilities  versus  range  are  displayed. 

The  results  are  summarized  in  Chapter  VII,  and  recom¬ 
mendations  for  future  work  are  given.  The  recommendations 


include  needed  developments  in  both  methodology  and  applica- 


tions.  The  greatest  need  for  applications  is  for  standar¬ 
dized  atatiatieal  descriptions  of  nnclear  effects  envi¬ 
ronment  s  . 

The  development  of  the  primary  numerical  method  nsed 
in  the  applications  problems  is  detailed  in  Appendix  A.  The 
documentation  of  several  modifications  to  a  nonparame tr ic 
estimation  method  is  included. 

A  source  code  listing  of  the  computer  program  NOSVET  is 
presented  in  Appendix  B. 

In  Appendix  C,  a  new  method  is  presented  for  finding 
the  distribution  of  monotonic  functions  of  random  variables. 
Functions  of  a  single  random  variable  are  considered  first, 
and  examples  presented.  The  theory  is  then  extended  to 
functions  of  two  random  variables,  and  tvo  examples  of  this 
are  given.  The  theory,  with  some  restrictions,  is  extended 
to  functions  monotonic  in  N  random  variables. 

In  Appendix  D,  the  numerical  technique  used  to  calcu¬ 
late  the  integral  of  conditional  density  functions  is 
explained  and  documented.  The  reliability  interference 
theory  integral  is  a  special  case  of  this  type  of  problem. 

Further  consideration  of  the  vulnerability  of  a  box- 
beam  wing  structure  in  nuclear  thermal  environments  is  pre¬ 
sented  in  Appendix  E.  In  ^articular,  the  relative  vulnera¬ 
bilities  of  the  parts  of  the  beam  are  discussed. 

Finally,  a  glossary  of  acronymns  and  nomenclature  is 
presented  in  Appendix  F,  keyed  by  Chapter. 
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Method! 


The  problem  to  be  considered  is  this:  An  aircraft  flys 
in  the  vicinity  of  a  unclear  weapon  at  the  time  that  it 
detonates.  If  the  aircraft  position  with  respect  to  the 
burst  is  known,  can  one  predict  the  probability  of  aircraft 
failure  due  to  the  effects  of  the  nuclear  weapon?  Vhat  is 
the  justification  for  such  a  prediction?  The  literature  is 
reviewed  below.  This  literature  search  leads  to  a  preferred 
approach  to  solving  this  problem. 


rminist ic  Methods 

The  Defense  Nuclear  Agency  (DNA)  document  DNA2048-H, 


Handbook  F 


on  Aircra 


[68]  is  representative  of  deterministic  methods  for 
assessing  aircraft  survivability.  The  basic  objective  of 
the  handbook  is  the  determination  of  sure-safe  (SS)  and 
sure-kill  ( SK )  ranges.  SS  and  SK  ranges  correspond  to 

range  locations  where  SS  and  SK  responses  of  an  aircraft 
subsystem  are  expected.  A  SS  response  corresponds  to 

incipient  damage,  or  limit  loads.  A  SK  response 
corresponds  to  a  catastrophic  damage  condition.  This  way 
of  thinking  has  its  origins  in  aircraft  design.  SS  loads 
are  limit  loads,  while  SK  loads  are  ultimate  loads,  or 
loads  beyond  ultimate.  Ultimate  loads  are  often  taken  to  be 


limit  loads  times  a  factor  of  safety,  typically  1.5.  This 
approach  to  survivability  allows  one  to  find  those  regions 
of  space  where  the  aircraft  is  snre-safe  or  snre-killed. 
Typical  outputs  of  the  method  of  Reference  [68]  are  shown  in 
Figures  2.1  and  2.2.  These  pictures  require  some 
expl anat ion. 

In  Figure  2.1,  the  target  occupies  the  origin  of 
coordinates.  A  nuclear  detonation  may  be  placed  anywhere  in 
the  plane  of  the  figure.  Any  burst  that  falls  on  the  solid 
contour,  or  outside  of  it,  leaves  the  target  at  the  origin 
completely  undamaged.  Hence,  the  region  of  space  outside  of 
the  contour  is  called  the  'sure-safe'  region.  The  probabil¬ 
ity  of  damage  is  nonzero  (but  unspecified)  in  the  region  of 
space  inside  the  contour. 

In  Figure  2.2,  the  target  also  occupies  the  origin. 
In  this  case,  however,  any  nuclear  burst  that  falls  on  the 
contour  or  inside  of  it  results  in  the  destruction  of  the 
target.  Hence,  the  region  inside  the  contour  is  now  the 

'sure-kill'  region.  Outside  of  the  contour,  the  probability 
of  damage  is  unspecified  (but  it  is  not  unity). 

The  chief  disadvantage  of  the  technique  is  that  it  is 
ambiguous.  There  are  regions  of  space  left  over  where  the 
state  of  the  aircraft  subsystem  cannot  be  determined. 
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The  obvious  problem  with  the  deterministic  technique 
can  be  eliminated  if  one  recognizes  that  the  failure  proba¬ 
bility  should  vary  smoothly  from  some  lov  value  in  the  SS 
region  to  some  high  value  in  the  SK  region.  A  number  of 
approaches  have  been  used  to  attempt  this,  and  they  are 
reviewed  below. 

A  i.  rc^r  a_f  t.  S.£  r  u£,£  u£  e„s  .  Bridgman  [44],  in  an  early 
Technical  Note,  attempted  to  find  a  continuous  probability 
of  damage  curve  by  postulating  a  damage  function  as  a  cumu¬ 
lative  lognormal  distribution  in  the  range  space.  This 
approach  closely  parallels  that  of  AP-5S0  [53],  discussed 
later.  Since  the  lognormal  is  a  two  parameter  distribution, 
an  infinite  number  of  possibilities  exists  for  the  contin¬ 
uous  damage  curve.  Bridgman  took  the  novel  approach  of 
coupling  the  information  in  DNA-2048H  so  as  to  find  the 
optimum  continuous  damage  model.  He  did  this  by  postulating 
the  damage  probability  to  be  some  arbitrarily  low  number 
(.02)  at  a  known  SS  range,  and  an  arbitrarily  high  number 
(.98)  at  a  known  SK  range.  The  SS  and  SK  ranges  were  deter¬ 
mined  by  the  methods  of  DNA-2048H.  Such  a  specification 
leads  to  a  unique  solution  for  the  two  parameters  of  the 
lognormal  distribution,  and  thus  specifies  the  damage  func¬ 
tion  . 

In  some  later  papers  [45,46],  Bridgman  recognized  that 
the  range  specifications  depended  not  just  on  the  target. 


bat  on  the  entire  soarce 


transmission,  and  interaction 


problem  as  veil.  He  then  postulated  a  continnons  failure 
distribution  in  the  stress  space,  again  finding  the  distri¬ 
bution  by  coupling  the  information  in  DNA-2048H  to  the 
lognormal  model  so  as  to  specify  the  survivability  function. 

There  have  been  other  attempts  to  treat  survivability 
as  a  statistical  problem  at  the  component  level.  The 
Studies  and  Analysis  group  at  Hq  DSAF  [42]  used  an  approach 
similar  to  Bridgman's.  Sure-safe  and  sure-kill  specifica¬ 
tions  in  the  stress  space  vere  assigned  as  percentiles  of 
the  component  failure  distribution.  The  possibilites  using 
parametric  statistics  vere  surveyed  as  shovn  in  Figure  2.3. 
Other  contributors  to  this  field  include  Gragg  [63].  The 
probabilistic  survivability  literature  for  aircraft  systems 
is  quite  limited  compared  to  that  of  ground  systems,  vhich 
vill  be  considered  next. 

Ground  Structures .  Ground  system  survivability  has 
been  studied  for  a  considerably  longer  time  than  has  air¬ 
craft  survivability.  One  of  the  earliest  probabil ist ic 
treatments  is  the  method  of  AP-550  [53].  The  failure  proba¬ 
bility  of  a  structure  is  given  as  a  complementary  cumulative 
lognormal  function,  vith  the  parameters  of  the  lognormal 
depending  on  the  target  hardness.  This  method  is  still  in 
use  by  the  targeting  community,  as  represented  by  the 
Defense  Intelligence  Agency  (DIA) .  In  this  method  one 
models  damage  directly  as  a  function  of  range,  based  on 


I 


Wt 


early  observations  of  tbe  Hiroshima  blast  damage.  Targets 
are  coded  by  a  special  number  as  being  either  overpressure 
sensitive  or  dynamic  pressure  sensitive.  The  special  coding 
of  a  target  type  yields  parameters  for  the  determination  of 
the  continuous  damage  function. 

A  quite  different  probabilistic  method  was  developed  by 
the  Defense  Nuclear  Agency  (DNA)  in  the  late  1960's  and 
applied  to  problems  of  Minuteman  Missile  survivability. 
This  technique  was  developed  independently  of  the  method 
used  by  DIA.  DNA  published  the  FAST  (Failure  Analysis  By 
Statistical  Techniques)  code  in  1974  [57].  FAST  is  an  ambi¬ 
tious  Monte  Carlo  code,  including  correlated  nuclear  weapon 
environments,  variable  system  reliability  models,  and  arbi¬ 
trarily  input  component  fragilities,  or  subsystem  failure 
distributions.  The  publication  of  FAST  was  a  significant 
milestone.  It  was  one  of  the  first,  if  not  the  first. 
Department  of  Defense  (DOD)  products  that  treated  surviva¬ 
bility  along  the  lines  of  classical  reliability  theory. 
After  being  used  in  the  Nuclear  Hardness  and  Evaluation 
Program  (NHEP)  [58],  it  was  not  used  much  until  1980.  At 
that  time  a  group  at  Lawrence  Livermore  National 
Laboratory  (LLNL)  began  using  it  to  investigate  probabilis¬ 
tic  modeling  of  survivability  to  EMP  [40,41]. 

Recent  probabilistic  survivability  studies  include 
examination  of  MX  survivability  in  both  its  horizontal  and 
vertical  shelters  [55].  These  studies,  also  conducted  by 


DNA,  did  not  use  FAST,  but  a  conbination  of  approximate 
probabilistic  techniques  [37]  and  Monte  Carlo  simulation 
[55].  Other  studies  of  hardening  of  ground  facilities  may 
be  found  in  the  literature  [49,50,75]. 

Another  notable  study  conducted  recently  by  DNA  is  that 
of  Jordano  [67].  He  investigated  the  uncertainty  of  nuclear 
effects  on  low  frequency  (LF)  communication  links.  He  cal¬ 
culated  the  probability  of  LF  signal  loss,  based  on  statis¬ 
tical  inputs  to  nuclear  debris  cloud  stabilization  altitude 
and  location.  This  work  differed  both  from  the  FAST  and  the 
approximate  probabilistic  techniques.  Rather  than  make 
distr ibut ional  assumptions,  Jordano  directly  calculated  the 
distribution  of  a  function  of  a  random  variable  based  on  a 
variable  transform  method. 

To  summarize,  ground  system  survivability  is  done  in 
two  ways.  The  targeting  community,  as  represented  by  D1A, 
treats  all  damage  as  an  empirical  function  of  range  [53]. 
The  nuclear  effects  community,  as  represented  by  DNA,  has 
developed  some  probabilistic  damage  models  based  on  classi¬ 
cal  reliability  theory  [37,55,57]  and  direct  variable  trans¬ 
formation  techniques  [67].  The  DNA  publications  just  men¬ 
tioned  form  the  nucleus  of  modern  probabilistic  survivabil¬ 
ity  literature  for  the  DOD.  However,  there  has  been  a  good 
deal  of  work  done  outside  of  DOD,  and  some  of  this  will  be 
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Other  Structure  t .  Two  groups  in  the  engineering  eosst- 
nity  outside  of  DOD  here  extensively  studied  aethods  for 
determining  the  probsbility  of  failure  of  structures.  Civil 
engineers  have  been  conducting  research  in  the  area  since 
the  late  1960's.  Most  of  the  literature  here  involves  the 
method  of  propagation  of  moments  of  random  variables,  abont 
which  more  will  be  said  in  Chapter  III.  This  technique  has 
been  applied  to  a  broad  range  of  structural  reliability 
problems  [7,27,28,30,31].  Some  of  these  even  deal  with 
nuclear  blast  damage  [71].  Textbooks  on  the  subject  of 
structural  reliability  are  beginning  to  be  published 
[15,23].  Since  about  the  early  to  mid  1970's,  the  nuclear 
engineering  community  has  also  taken  an  active  interest  in 
the  subject.  Proceedings  of  the  conferences  on  Structural 
Mechanics  in  Reactor  Technology  (SMIRT)  include  many  papers 
on  the  application  of  probabilistic  methods  to  problems  of 
nuclear  safety.  Here  too,  books  are  being  written  [25,29]. 
In  most  cases,  these  methods  and  applications  are  based  on 
mathematical  reliability  theory,  coupled  with  standard  prop¬ 
agation  of  moments  techniques. 

Summary 

The  literature  just  reviewed  includes  deterministic  and 
probabilistic  methods.  Deterministic  methods  do  not  provide 
for  continuous  failure  (or  survivability)  functions.  This 
violates  intuition.  On  the  other  hand,  probabilistic  models 
that  are  based  on  direct  modeling  with  range  are  justified 
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Overview 


The  mathematics  of  reliability  theory,  and  the  field  of 
reliability  analysis,  came  about  as  a  result  of  problems  the 
DOD  had  with  electronics  piece-parts  in  the  early  1940's 
[8].  Only  since  about  the  mid  1960's  has  the  field  been 
extended  to  problems  of  large-scale  engineering  systems.  In 
the  paragraphs  that  follow,  the  basic  theory  of  reliability 
mathematics  is  reviewed.  Engineering  determinism  is  shown 
to  fit  naturally  into  the  theory  as  a  special  case.  The 
problems  in  applying  the  theory  to  engineering  systems  are 
reviewed,  as  are  previous  approaches  to  solving  those 
problems . 

Mathematical  Reliability  Theory 

Strength  Distributions.  The  strength  distribution  is  a 
statistical  property  of  a  population  of  devices  operating  in 
some  stress  environment.  Examples  include  light  bulbs 
operating  for  a  given  number  of  hours,  or  samples  of  wing 
materials  under  tensile  testing.  The  strength  distribution 
is  sometimes  called  the  failure  distribution  or  the 
resistance  distribution.  It  is  a  measure  of  how  well  the 
devices  resist  a  stress  load.  The  distribution  is  most 
easily  determined  for  the  case  when  a  number  of  items  can 
all  be  tested  to  failure.  In  that  case,  one  can  plot  the 
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fractional  number  of  failures  occuring  in  a  given  stress 
interval  versus  the  midpoint  of  the  interval.  The  resulting 
histogram  approximates  the  underlying  probability  density 
function  (PDF)  for  failure.  This  statistical  information  is 
a  measure  of  the  strength  of  the  item,  since  the  central 
value  of  the  distribution  indicates  the  value  of  stress 
where  most  failures  occur.  In  addition,  the  amount  of 
scatter  about  the  central  value  measures  the  variability  in 
the  quality  of  the  item.  The  stress  value  at  which  an  item 
will  fail  is  a  random  variable.  The  probability  that  an 
item  will  fail  at  a  particular  stress  value  s  is  then 
expressed  by: 

p  f  =  Pr{Sls}=y^oof  s(s)ds=Fg(s)  (3.1) 

where  pf  is  the  failure  probability,  S  the  strength  random 
variable,  s  a  particular  value  of  S,  fg(s)  the  strength  PDF, 
and  Fg(s)  the  strength  cumulative  distribution  function 
(CDF). 

If  the  testing  of  n  items  resulted  in  every  single  item 
failing  precisely  at  the  value  SK,  then  a  true  cookie-cutter 
distribution  would  result.  The  PDF  would  be  a  Dirac  delta 
function  [2]  located  at  SK,  and  the  CDF  would  be  a  step 
function  as  shown  in  Figure  3.1  (top).  A  more  likely  result 
is  that  every  item  fails  at  a  slightly  different  value, 
resulting  in  a  CDF  similar  to  that  of  Figure  3.1  (bottom). 


Cookie-Cutter  Distribution 


Uniform  Distribution 


COMPONENT  RESPONSE 


Figure  3.1.  Cookie-Cutter  (Top)  and  Random  (Bottom) 
Failure  Distributions 

(Adapted  From  Reference  [57]) 


Si£e.s.s.  Pis  tribat  ions.  Equation  3.1  gives  tbe  failure 
probability  for  the  item  provided  the  stress  is  exactly 
known.  As  in  the  case  of  the  failure  point  of  the  device, 
the  applied  stress  may  not  be  known  exactly.  Repeated 
measurements  of  the  free  field  stress  may  yield  a  central  or 
expected  value  f  j  r  the  measurements,  but  also  dispersion 
about  the  central  value.  The  applied  stress,  s,  would  then 
also  be  a  random  variable,  described  by  a  statistical  dis¬ 
tribution  . 

The  Interference  Integral  and  Related  Random  Variables. 
If  the  distributions  of  S  and  a  have  been  determined,  then 
the  general  expression  for  the  failure  probability  is  given 
by : 

Pf=Pr(Sls}  (3.2) 

Provided  the  strength  and  stress  variables  are 
independent.  Equation  3.2  can  be  written  as: 

Pr(S<s}=/®  fs(s)Fs(s)ds  (3.3) 

QD 

where  p£  is  the  failure  probability  as  before,  S  the 
strength  random  variable,  s  the  stress  random  variable, 
fs(s)  the  PDF  of  the  stress  distribution,  and  Fg(s)  the  CDF 
of  the  strength  distribution.  In  words,  the  failure 
probability  is  just  the  chance  that  a  random  selection  from 
the  sample  space  of  all  strength  values  yields  a  number  less 
than  or  equal  to  a  random  selection  from  the  sample  space  of 
all  stress  values.  Equation  3.3  is  often  referred  to  as  the 
reliability  interference  integral,  or  just  the  interference 
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integral,  since  contributions  to  the  integral  occur  only  in 


the  variable  space  where  the  density  functions  overlap  or 
“interfere”.  The  shaded  portion  of  Figure  3.2  illustrates 
the  interference  region.  Equation  3.3  is  seen  to  reduce  to 
Equation  3.1  if  the  stress  PDF  is  a  Dirac  delta  function. 

If  the  stress  and  strength  random  variables  are  not 
independent.  Equation  3.3  cannot  be  used  as  written. 
However,  Equation  3.2  is  still  valid,  and  the  failure 
probability  can  be  calculated  by  finding  the  distribution  of 
either  the  margin  variable,  (  ,  or  the  safety  factor 
variable,  17  .  The  margin  variable  is  defined  by: 

£  =  S-s  (3.4) 

» 

Consequently,  in  terms  of  the  margin  variable,  the  failure 
probability  could  also  be  calculated  by: 

Pf  =  Pr{$10}  (3.5) 

Also,  the  safety  factor  variable,  provided  the  applied 
stress  is  not  zero,  is  defined  by: 

n  =  S/s  (3.6) 

In  terms  of  the  safety  factor,  the  failure  probability  can 
be  calculated  by: 

Pf  =  Prl»?i.l}  (3.7) 

At  this  point,  examination  of  Equations  3.4  and  3.6 
shows  that  one  of  the  requirements  for  solving  the  problem 
is  the  ability  to  find  the  distribution  of  a  function  of  two 
random  variables.  This  can  be  a  difficult  problem,  and  it 
will  be  discussed  in  more  detail  below.  Before  considering 
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that,  a  discussion  of  the  reduction  of  the  interference 
integral  for  the  deterministic  case  is  presented. 

Reduction  To  Engineering  Determinism.  If  the  stress  is 
deterministic,  the  PDF  of  a  may  be  represented  by  a  Dirac 
delta  function  centered  at  the  known  stress  value,  say  "s. 
That  is, 

fs(s)=$(s-?)  (3.8) 

Using  Equation  3.8  in  the  interference  integral  of  Equation 
3.3  reduces  the  integral  to: 

pf  =  Fs(i)  (3.9) 

If  the  strength  is  also  deterministic,  then  the  strength  CDF 
is  given  by: 


Fs(s)=0  v  s  <  SK  (3.10) 

Fs(s)=l  r  s2SK  (3.11) 

The  results  of  Equations  3.10  and  3.11  are  the  familiar 
deterministic  cookie-cutter  damage  functions. 

Going  the  other  way.  Equations  3.10  and  3.11  can  be 
used  first  inside  the  interference  integral  of  Equation  3.3. 
Since  the  integrand  is  zero  for  every  s<SK,  the  limits  are 
changed  so  that: 

(3.12) 

But  Equation  3.12  is  just  the  statement  that: 
p  f  =  Pr  (a  >SK)  (3.13) 

Again,  by  hypothesis,  only  one  value  of  a  is  ever 
observed,  namely  s.  Consequently,  Equation  3.12  takes  only 


Pf=  f,(s)ds 

•'sk 


two  values: 


p£*l  ▼  s2.SK 
Pi=0  ▼  s<SK 


(3.14) 

(3.15) 


The  set  of  Equations  3.14  and  3.15  is  identical  in 
meaning  to  the  set  of  Equations  3.10  and  3.11.  Cookie- 
cotter  damage  distributions  are  the  result. 

The  reduction  of  the  general  reliability  integral 
approach  to  the  c ooki e- c u t t e r  case  for  Dirac  delta  PDF's 
makes  it  a  natural  choice  for  adding  probabilistic 
information  to  the  survivability  problem.  One  sees  that 
engineering  determinism  is  nothing  more  than  a  subset  of  a 
more  general  probabilistic  approach.  The  subset  is  found  by 
replacing  all  continuous  probability  density  functions  by  a 
very  special  probability  density  function — the  Dirac  delta, 
which  has  a  single  parameter.  Seen  in  this  light, 
engineering  determinism  would  seem  to  have  no  inherent 
theoretical  superiority  over  probabilistic  methods  based  on 
reliability  theory.  Even  so,  probabilistic  approaches  to 
survivability  problems  are  not  common.  The  reasons  for  this 
are  examined  next. 


Applications  Difficulties  For  Engineering  Systems 

The  direct  application  of  s tr e s s- s tr e ng th  interference 
theory  has  been  thwarted  by  problems  in  system  reliability 
modeling,  analytic  and  numerical  intractability  in  finding 
the  distribution  of  functions  of  random  variables,  and  a 
serious  lack  of  hard  data.  These  are  all  discussed  in  turn 
below. 
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Sx.£t.e.m  Ke._lia.bjL_lj._tx  M.°  4«.  i-lfi-g. •  Even  if  data  were 
abundant  and  mathemat  leal  methods  available,  it  is  still 
difficult  to  determine  bow  individual  component  failures  can 
combine  to  create  a  system  failure.  An  example  from  the 
aircraft  survivabil ity  problem  is  the  melting  or  rupturing 
of  non- s  t  r  u  c  t  ur  a  1  skin  panels.  It  simply  is  not  true  that 
the  loss  of  a  single  skin  panel  will  mean  the  loss  of  the 
aircraft.  It  is  true  that  as  panels  are  lost,  the  drag 
coefficient  increases,  and  that  at  some  point,  as  more  and 
more  panels  are  lost,  the  drag  coefficient  is  such  that  the 
mission  cannot  be  completed — i.e,  the  system  has  failed. 
The  difficulty  in  constructing  a  model  of  system  reliability 
or  survivability  as  a  function  of  component  survivabil ity  is 
apparent.  System  reliability  modeling  is  discussed  in  most 
reliability  engineering  textbooks  [8]. 

An  approach  used  in  reactor  safety  studies  is  that  of 
fault  tree  analysis  [29].  Basically,  a  list  is  created  of 
every  event  that  is  to  be  the  subject  of  analysis.  These 
events  are  termed  top  events.  A  functional  diagram  of  the 
system  is  then  studied  so  as  to  identify  contributing  events 
that  may  directly  cause  the  top  event  to  occur.  These 
contributing  events  typically  cause  the  top  event  to  occur 
through  AND/OK  Boolean  operations.  For  example  [29],  a 
circuit  breaker  may  fail  to  trip  because  it  fails  randomly 
by  itself,  OR  a  trip  signal  is  not  received.  The  con¬ 
tributing  events  can  also  be  further  analyzed.  To  continue 


the  example,  perhaps  the  trip  signal  is  not  received  because 


relay  1  AND  relay  2  remained  closed.  Fanlt  tree 
constrnction  can  become  a  substantial  task  for  a  large 
engineering  system. 

At  this  point,  it  should  be  noted  that  whatever  number 
is  obtained  by  a  probabilistic  assessment  of  survivability, 
it  will  not  represent  an  absolute  measure  of  the  frequency 
of  system  failure.  Since  no  system  reliability  diagram  can 
be  perfect,  a  probabilistic  estimate  of  survivability  is  a 
conditional  probability  measure.  At  the  very  minimum,  it  is 
conditional  due  to  the  system  reliability  model.  Since  it 
is  not  possible  to  know  all  the  failure  mechanisms  of  a 
system,  the  resulting  survivability  number  is  an  upper  bound 
to  the  system  survivability.  Thus,  for  military 
applications,  targeters  have  more  to  gain  by  this  type  of 
analysis  than  do  defenders. 

Mathemat ical  Methods.  Another  difficulty  in  the  direct 
application  of  reliability  theory  is  the  analytic  and 
numerical  intractability  of  the  mathematics  of  statistics. 
The  most  common  problem,  and  one  already  hinted  at,  is  that 
of  finding  the  distribution  of  a  function  of  one  or  more 
random  variables.  Very  few  problems  have  exact  analytic  or 
n  1  AtAd-fnrm  tAfntinnc.  IIa  thodi  t  hi  t  hi  v  a  Ha  in  need  b  v 


survivability  analysts  include  propagation  of  moments  with  a 
priori  distributional  assumptions,  Monte  Carlo  simulation, 
and  variable  transformation  techniques. 
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Propagation  of  moments  methods  have  been  nsed  by  DNA  in 
some  of  the  studies  of  missile  basing  [55].  Basically,  given 
an  equation  of  the  form: 

Z=g(I1.X2.  .  .  .XB)  (3.16) 

moment  information  about  the  X^  is  collected  (usually  the 
mean  and  variance).  The  mean  and  variance  of  Z  are  then 
calculated  by  means  of  Taylor  Series  expansions  about  some 
point  [8,12].  This  requires  the  determination  and 
evaluation  of  partial  derivatives.  Once  the  mean  and 
variance  of  Z  have  been  estimated,  a  lognormal  distribution 
function  for  Z  is  usually  assumed  [71].  Correlations  among 
the  Xj  must  be  accounted  for.  Along  these  same  lines,  a  new 
technique  based  on  a  linear  algebra  has  been  developed  by 
Ditlevsen  [23].  This  approach  assumes  that  only  moment 
information  is  available  for  all  inputs,  and  thus  only 
moment  information  can  be  used  for  all  outputs.  Ditlevsen's 
linear  algebra  approach  seeks  to  accomplish  the  same  goals 
as  a  probability  theory,  yet  can  be  presented  without  a 
reference  to  probability  theory. 

There  are  several  disadvantages  to  using  propagation  of 
moments.  One  is  that  not  all  combinations  of  random 
variables  yield  distributions  that  have  moments.  The 
Cauchy,  for  example,  cannot  be  estimated  efficiently  by  such 
methods  [23,13].  Another  disadvantage  is  the  common  use  of 
a  priori  distributional  assumptions.  However,  perhaps  the 
worst  disadvantage  is  the  tediousness  of  calculating  many 


I 


partial  derivatives.  For  example,  in  Appendix  B  of 
Reference  [12],  one  can  find  an  expression  for  the  expected 
value  of  a  fourth  moment.  This  single  expression  occupies 
48  lines  of  typewritten  text!  Needless  to  say,  tasks  of 
this  magnitude  axe  prone  to  errors.  Most  of  these 
objections  can  be  overcome  by  using  Monte  Carlo  simulation. 

Direct  Monte  Carlo  simulation  is  described  by 
Stancampiano  [32].  This  is  basically  a  brute  force 
counting  operation.  One  samples  randomly  the  distribution 
of  margin  or  safety  factor.  This  is  done  by  randomly 
drawing  from  the  strength  and  stress  distributions,  which 
may  be  functions  of  multiple  random  variables,  and  counting 
the  relative  number  of  times  the  margin  variable  is  less 
than  0,  or  the  safety  factor  less  than  1.  If  the  part  has  a 
very  high  reliability,  this  counting  operation  can  become 
expensive  . 

Monte  Carlo  with  fitting  seeks  to  overcome  the  expense 
of  brute  force  Monte  Carlo  counting.  In  this  case,  one  fits 
the  simulated  data  to  a  distribution  function,  perhaps  a 
normal,  or  some  other  type.  Stancampiano  shows  [32]  that 
some  problems  can  be  sensitive  to  the  choice  of  distribution 
function  chosen.  An  elegant  technique  is  to  calculate  high 
order  moments  of  the  data  and  use  the  Shannon  maximum 
entropy  method  to  find  the  minimally  prejudiced  probability 
distribution  [13].  This  is  not  always  straight-forward, 
since  the  solution  for  the  minimally  prejudiced  PDF  involves 
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solving  a  non-linear  programming  problem.  Also,  there  is  no 
a  priori  guarantee  that  the  distribution  of  the  output 
variable  has  moments. 


Other  Monte  Carlo  techniques  exist,  though  they  will 


Importance  sampling  [32]  can  be  used  if  one  knows  that  the 
location  parameter  for  the  stress  distribution  is  relatively 


efficiency  by  concentrating  the  sampling  in  the  area  of 
interference  between  the  strength  and  stress  density 
functions.  This  will  not  be  useful  for  the  problem  at  hand 
since  the  applied  stress  from  a  nuclear  weapon  cannot,  in 
general,  be  confined  to  a  limited  range  of  values. 

Other  methods  include  the  use  of  specialized 
mathematics,  such  as  Cook's  work  with  B  Functions  [4]. 
These  methods  are  still  cumbersome  and  unattractive  at 
present  for  engineering  applications. 

S  pii  p  s.  e,  D a jt a  Se.  _£  s. .  Finally,  the  lack  of  data  for  the 
required  inputs  to  a  probabilistic  assessment  limits  the 
usefulness  of  the  reliability  theory  approach.  The  most 
serious  limitation  of  the  lack  of  data  is  in  the  ability  to 
identify  parametrically  the  input  distributions.  With  a 
small  sample,  it  is  possible  that  a  number  of  density 
functions  would  slip  through  even  the  best  of  goodness  of 
fit  tests.  Ashley  [3]  has  shown  how  the  use  of  a  weak  test, 
such  as  the  Chi-square  which  is  popular  among  engineers,  can 


how  many  there  are. 

The  dearth  of  hard  data  has  led  aan;  to  consider 
Bayesian  estimation  techniques.  The  use  of  Bayesian 
statistics  in  reliability  theory  is  discussed  by  Kapur  and 
Lamberson  [8].  Bayesian  statistics  combines  subjective 
judgment  or  experience  with  hard  data  to  provide  atatistical 
information  similar  to  the  classical  statistical  inference 
approach.  The  method  is  controversial  since  “no 
experimental  or  analytical  methods  exist  for  the 
quantification  of  .  .  .belief"  [8].  Nevertheless,  Bayesian 
methods  have  been  widely  used  in  doing  probabilistic  risk 
assessments.  The  reasons  for  this  are  best  represented  by 
direct  quotes  from  the  literature: 

Ang  writes  in  1975,  "When  the  observed  data  are 
limited,  as  is  often  the  case  in  engineering,  the 
statistical  estimates  have  to  be  supplemented  (or  even 
superseded)  by  judgmental  information.  With  the  classical 
statistical  approach,  there  is  no  provision  for  combining 
judgmental  information  with  observational  data  in  the 
estimation  of  parameters."  [1] 
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McCormick,  in  a  1981  book  on  risk  analysis  [29],  writes 
"The  interpretation  of"  (the  law  of  large  numbers)  "is  clear 
enough  for  experiments  that  can  be  repeated.  There  are  many 
occasions  in  which  the  knowledge  available  is  less  precise, 
especially  when  the  engineer  deals  with  rarely  occurring 
events  that  form  the  basis  of  many  risk  evaluations.  Then 
it  is  necessary  to  resort  to  the  axiomatic  or  subjective 
approach  to  the  concept  of  probability,  which  we  shall  use 
from  now  on.  " 

Subjective  probabilistic  assessment  is  also  mentioned 
by  Bevensee  [41].  The  lack  of  hard  data  for  EMP 
survivability  assessment  is  apparently  so  serious,  that  in  a 
1981  follow-up  study  to  that  of  Reference  [41],  fully  29 
pages  of  a  draft  Lawrence  Livermore  National  Laboratory 
(LLNL)  report  are  devoted  to  methods  of  polling  expert 
opinion  [39].  Personal  conversations  with  experts  in  the 
study  of  ground  systems  survivability  confirm  that  opinion 
surveys  are  a  common  source  of  information  [48]. 

This  lack  of  hard  data,  combined  with  the  controversy 
surrounding  the  use  of  Bayesian  statistics,  leads  to  a 
perplexing  dilemma.  Faced  with  the  difficulties  in 
rigorously  pursuing  an  approach  along  classical  statistical 
lines,  one  may  wish  to  retain  the  traditional  approach  of 
engineering  determinism.  On  the  other  hand,  engineering 
determinism  cannot  deal  with  uncertainty  other  than  by  worst 
case  analysis  or  by  parametric  surveys.  From  a  probability 
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theory  point  of  view,  one  has  not  avoided  subjectivity  by 
retaining  engineering  determinism.  One  has  simply  chosen 
the  Dirac  delta  function  as  the  underlying  PDF  for  all 
variables  of  interest.  The  present  probabilistic  alter¬ 
natives  to  that  are  to  choose  continuous  distributions,  and 
somehow  to  estimate  the  parameters  for  those  distributions. 

Summary 

To  summarize  the  Chapter  then,  three  of  the  more 
serious  difficulties  in  applying  stress-strength 
interference  theory  to  engineering  systems  have  been  noted. 
These  include:  (a)  the  problem  of  system  survivability 
modeling  based  on  component  survivability,  (b)  the 
difficulty  inherent  in  the  mathematical  methods  for 
propagating  statistical  information,  and  (c)  the  sparseness 
of  the  databases  for  engineering  inputs. 

In  Chapter  IV  (with  related  appendices),  a  new  approach 
for  dealing  with  problems  (b)  and  (c)  above  will  be 
presented.  It  will  be  shown  that  newly  developed  methods  in 
nonp a r a m e t r i c  estimation  provide  a  tool  for  finding  the 
distribution  of  a  function  of  multiple  random  variables.  It 
will  also  be  argued  that  nonparametr ic  estimation  provides  a 
logical  method  for  propagating  uncertainty  in  an  engineering 
mode  1 . 
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Tool 


IV.  Noacttiiottlc  Estl»ation-A  New  |atg££tam  IhSSLSJL 

Overview 

In  this  Chapter,  an  effort  will  be  made  to  deal  with 
some  of  the  problems  mentioned  in  Chapter  III.  However,  the 
problem  of  system  survivability  modeling  from  knowledge  of 
component  survivability  will  not  be  dealt  with.  In  the 
chapters  that  follow,  an  aircraft  system  will  be  treated  as 
a  simple  series  reliable  assembly.  However,  much  more  can 
be  done  in  overcoming  the  problems  involved  with  mathemati¬ 
cal  methods,  and  with  sparse  data  bases.  These  are  pre¬ 
sented  at  this  time. 

Distribution  o  f  a.  Funct  ion  of  Multiple  Random  Variables 

Some  recent  developments  [14]  in  nonparametric  estima¬ 
tion  techniques  provide  a  new  and  powerful  statistical  tool 
for  use  in  survivability  calculations.  This  tool  can  be 
used  as  an  empirical  one,  for  use  on  observed  sets  of  data. 
It  can  also  be  used  as  a  numerical  method  in  its  own  right. 
Although  most  of  the  details  are  discussed  in  Appendixes  A 
and  C,  a  basic  outline  of  the  approach  is  discussed  below, 
and  an  example  from  the  reactor  risk  analysis  literature  is 
presented. 

As  noted  earlier,  one  of  the  tasks  that  must  be  done  in 
order  to  make  s t r e s s- s t r e n g t h  interference  theory  work  is 


that  of  finding  the  distribution  of  a  function  of  one  or 


more  random  variables.  With  the  nonparametr ic  estimation 
technique,  the  following  can  be  demonstrated: 

(a)  Given  the  random  variable  equation  Z=g(X),  with 
g(x)  monotonic,  and  the  distribution  function  of  X  (whether 
parametric  or  non— parametric),  the  distribution  function  of 
Z  can  be  exactly  determined  at  selected  points. 

(b)  Given  the  random  variable  equation  ZBg(X,X),  the 
distribution  functions  of  X  and  X,  and  the  restriction  that 

vx,  dg/^y/O  vy,  then  the  distribution  function  of  Z 
can  be  determined  at  a  selected  number  of  points  by  using 
the  equation  of  conditional  probability  as  a  sampling  rule. 
The  accuracy  in  the  distribution  function  of  Z  is  limited 
only  by  standard  problems  of  numerical  integration. 

(c)  Given  the  random  variable  equation 
Z=g(X1.X2.  .  .  .  Xa)  (4.1) 

the  distribution  functions  of  the  X^,  and  the  restriction 
that  dg/dx^fO  v  x ^ ,  then  the  distribution  function  of  Z  can 
be  determined  at  a  selected  number  of  points,  provided  the 
X^  are  independent  or  have  known  correlations.  The  accuracy 
in  the  resulting  distribution  function  of  Z  is  limited  only 
by  standard  problems  of  numerical  integration. 

This  technique  is  best  illustrated  by  working  an 
example  from  the  actual  survivability  literature.  The 
density  function  for  the  safety  factor  of  a  reactor  pressure 
vessel  will  be  calculated  by  the  above  method,  and  compared 
to  a  standard  Monte  Carlo  calculation  [32]. 


A  Benchmark  Problem- Safety  Factor  o f  a  Reactor 
Pressure  Vessel 

A  very  long  ductile  steel  cylinder  is  subjected  to  an 
applied  internal  pressure  P.  The  applied  pressure  P  is 
found  to  be  a  random  variable  with  a  three-parameter  Weibull 
distribution.  That  is: 

Fp (p ) =1- e xp  [- ( (p- c ) / a)b ]  y  p>c  (4.2) 

Fp(p)=0  ▼  p^c  (4.3) 

where  p  is  a  particular  pressure  value  in  EPa  ( k i 1  op a s c a  1 s ) , 
a  is  the  Weibull  scale  parameter  in  EPa,  b  the  dimensionless 
Weibull  shape  parameter,  and  c  the  Weibull  location  para¬ 
meter  in  EPa.  Stancampiano  gives  the  values  as: 


a=.10665P0 

(4.4) 

b  =  4 .212 

(4.5) 

c=.9P0 

(4.6) 

Pq=43 1 6  EPa 

(6  26  p  s i) 

(4.7) 

The  random 

variable  P  is  the  stress  variable. 

and  it 

has  the  parametric  representation  just  noted.  However,  one 
may  represent  P  nonpar am e t r i c a  1 ly  rather  than  parametric¬ 
ally.  This  is  done  by  solving  Equation  4.2  for  the  m  values 
of  p^  that  satisfy  the  equation: 

Pi  =  Fp1(Gi)  (4.8) 

where  the  are  given  by: 

Gi=( i+a) / (m+p)  ;  i  =  l,m  (4.9) 

and  a  and  P  are  constants  satisfying: 


-li*iPil 


(4.10) 


The  m  values  of  p  found  by  solving  Equation  4.8  are 
ordered  from  smallest  to  largest  and  denoted  by  the  set  of 
values: 

{ P i } i i=l 

As  discussed  in  Appendix  C,  if,  for  every  i  from  1  to 
m,  the  of  Equation  4.9  are  plotted  versus  the  ordered  p^ 
in  the  set  just  described,  then  the  result  in  general  is  an 
approximation  to  the  distribution  function  Fp(p).  However, 
for  this  case,  the  result  is  exact  at  the  data  points  p^. 
This  is  because  the  data  (the  p^)  were  drawn  by  the  rule  of 
Equation  4.8.  Such  a  set  of  points  will  be  defined  as  a 
stylized  set  (Reference  Appendix  C). 

The  distribution  of  P  derived  from  a  set  of  SO  stylized 
points  is  shown  in  Figure  4.1.  The  plotting  rule  used  was 
Hazen's  Rule  where  a  and  0  are  given  by  -.5  and  0  respec¬ 
tively.  The  distribution  function  of  Figure  4.1  is  exact  at 
the  percentiles  .01(.02).99  and  in-between  values  are  simply 
linearly  interpolated  from  these.  As  discussed  in  Appendix 
A,  the  end-point  values  are  found  by  requiring  the  integral 
of  the  PDF  to  be  unity. 

At  this  point,  it  might  appear  that  this  representation 
of  the  distribution  of  P  is  needlessly  complicated,  since 
(a)  the  distribution  of  P  is  already  known  parametrically, 
and  (b)  the  parametric  form  was  used  anyway  to  provide  the 
data  to  the  nonparametric  technique. 
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The  answer  to  the  first  objection  is  that,  although  it 
seems  complicated  now,  the  application  of  this  method  will 
make  more  difficult  problems  much  simpler  later  on.  To  find 
a  distribution  function  parametrically,  one  must  first  know 
the  density  function  or  distribution  function,  and  then  all 
the  parameters  that  go  into  it.  If  these  were  always  known 
for  every  conceivable  combination  of  random  variables,  then 
density  function  identification  would  be  greatly  simplified. 
As  it  is,  some  of  the  simplist  combinations  of  random  vari¬ 
ables  (Reference  Appendix  C)  have  complicated  density  and 
distribution  functions. 

As  to  the  second  objection,  the  fact  that  the  nonpara- 
metric  technique  is  using  data  drawn  from  a  parametric  model 
is  just  a  property  of  the  example.  Stancampiano  [32]  does 
not  say  how  the  distribution  of  P  was  found  to  be  three- 
parameter  Veibull.  If  one  simply  believes  him,  then  the 
parametric  input  is  accepted.  If  one  does  not  believe  him, 
then  one  obtains  the  data  for  himself,  and  estimates  the 
distribution  of  P  independently,  using  another  parametric 
model,  or  perhaps  a  no  np  a  r  a  m  e  t  r  i  c  one.  In  any  case,  every 
input  variable  will  have  its  distribution  function  specified 
somehow.  No  matter  how  it  is  done,  the  nonparame tr ic  repre¬ 
sentation  can  be  used. 

To  continue  the  calculation  of  the  survivability  of  the 
pressure  vessel,  one  now  needs  to  find  the  distribution  of 
the  strength,  having  determined  the  distribution  of  the 


applied  stress.  At  this  point,  a  subtle  departure  from  the 
norm  of  classical  reliability  theory  oecurs--a  departure 
shared  by  most  problems  in  nuclear  survivability.  Classic¬ 
ally,  the  strength  distribution  would  be  found  by  taking  a 
large  number  of  cylinders,  and  testing  them  all  to  failure. 
The  strength  distribution  would  then  be  determined  by  the 
standard  methods  of  statistical  inference.  However,  in  this 
case,  as  in  most  cases  involving  nuclear  survivability,  the 
items  in  question  are  far  too  few  in  number,  and  too  expen¬ 
sive  per  copy,  to  test  to  destruction.  Instead,  one  models 
the  response  of  the  part.  For  the  cylinder  example,  the 
point  at  which  the  vessel  is  believed  to  burst  is  given  by 
the  equation: 

Pb=Atf;FcyllnW  (4*11) 

where  A  is  a  dimensionless  modeling  coefficient,  <T^  is  the 

engineering  ultimate  tensile  strength  in  KPa,  Fcyl  is  a 

dimensionless  function  of  strain  and  elongation  given  by: 

Fcyl=,25/(<u+,227)(e/fu) €®  (4.12) 

is  the  true  strain  at  maximum  tensile  test,  and  is 
given  as  a  function  of  the  uniform  elongation  ,  ( 'n  by: 

tu=ln(l+c;)  (4.13) 

W=D/d  (4.14) 

d=,D-2 1  (4.15) 

t  being  the  wall  thickness  of  the  cylinder,  D  the  outer 
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diameter#  and  d  the  inner  diameter 


The  distribution  of  the  strength  variable,  that  is,  the 
distribution  of  P^,  is  determined  by  finding  the  distribu¬ 
tion  of  the  individual  input  variables,  and  then  finding  the 
distribution  of  the  funotioa  of  the  inputs.  The  input 
distributions  are  given  by: 

A  is  normally  distributed  with  mean  1.0478  and  standard 
deviation  .0948.  <7^  is  normal  with  a  mean  of  433,000  KPa 
and  standard  deviation  of  21,200  KPa,  is  normal  with  a 
mean  of  .4485  and  standard  deviation  of  .0377.  The  corre¬ 
lation  between  the  tensile  strength  and  elongation  is  an 
additional  complication.  The  correlation  coefficient  was 
found  to  be  -.498.  The  outer  diameter  D  was  taken  as  uni¬ 
formly  distributed  between  60.639  cm  and  61.281  cm.  The 
wall  thickness  t  was  taken  as  uniformly  distributed  between 
1.237  cm  and  1.532  cm  [32].  The  nonp  a  r  am  e  t  r  i  c  representa¬ 
tions  of  the  input  variables  are  shown  in  Figures  4.2 
through  4.6. 

The  strength  distribution  can  now  be  fonnd  by  a 
sequence  of  binary  operations.  Some  variable  changes  sim¬ 
plify  the  equation  for  the  burst  pressure,  P^.  The  geomet¬ 
rical  factors  can  be  represented  by: 

X1  =  ln¥  =  ln(D/(D-2t))  (4.16) 

The  distribution  of  X^,  found  by  the  method  described 
earlier,  is  shown  in  Figure  4.7.  With  the  distribution  of 
Xj  determined,  the  dimension  of  the  random  variable  equation 
is  reduced,  and  the  burst  pressure  can  be  written  as: 
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(4.17) 


To  handle  the  correlated  variables,  the  variable  X2  it 
defined  by: 

The  distribution  of  X2  is  shown  in  Figure  4.8.  At  this 
point,  the  random  variable  equation  for  the  burst  pressure 
has  been  reduced  to  an  independent  set  of  variables: 

Pb=AX2Xl  (4*19) 

One  can  now  proceed  with  a  sequence  of  binary  opera¬ 
tions  to  find  the  distribution  of  the  burst  pressure. 
Defining  the  variable  X3  by: 

X3=AX1  (4.20) 

the  equation  for  the  burst  pressure  is  reduced  to: 

Pb=x2x3  (4.21) 

The  distribution  of  X3  is  illustrated  in  Figure  4.9. 

At  this  point,  there  are  at  least  three  ways  to  pro¬ 
ceed.  One  could  find  the  distribution  of  Pb  from  Equation 
4.21  above,  and  then  directly  integrate  to  find  the  failure 
probability.  The  integration  technique  is  discussed  in 
Appendix  D.  It  is  analagous  to  the  graphical  estimation 
technique  [8],  or  the  Mellin  transform  technique  [9]. 
Alternately,  we  could  evaluate  the  distribution  of  the  mar¬ 
gin  variable,  £,  or  the  safety  factor  variable  T)  .  The 

latter  will  be  done.  In  this  case,  *7  is  given  by: 

T?-Pb/P-X2X3/P  (4.22) 
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Defining  X4  by: 

X4=X3/P  (4.23) 
the  distribution  of  X4  is  determined  and  illustrated  in 
Figure  4.10.  The  final  distribution  of  the  safety  factor, 
given  by: 

t)=X2X4  (4.24) 
is  illustrated  in  Figure  4.11. 

One  of  S t a nc a mp i a no ' s  single  random  samples  of  the 
safety  factor  is  displayed  in  Figure  4.12.  It  is  of  size 
3500.  In  Figure  4.13  the  results  of  the  nonparame tr ic 
technique  are  overlayed  with  that  of  the  random  sample 
method.  As  mentioned  earlier,  the  nonparametric  results  are 
exact  (within  integration  errors)  at  the  percentiles 

t 

.01 (.02)  .99.  Also,  the  nonparametric  technique  is  efficient 
enough  to  run  on  a  Z-80  nicroconputer.  The  results  of 
Figures  4.1  to  4.11  were  done  on  such  a  computer.  A  check 
was  made  on  the  results  by  doing  the  problem  on  a  CYBER  750 
also.  For  this  check,  100  stylized  points  were  used  to 
represent  each  distribution.  Except  for  increased  tail 
information,  the  results  were  identical  to  that  for  50 
stylized  points. 

Of  course,  one  need  not  use  the  nonparametric  algorithm 
in  such  a  rigorous  fashion  as  just  described.  Instead,  one 
may  use  random  Monte  Carlo  sampling  just  as  Stancampiano 
did,  and  then  use  the  nonparametric  method  on  the  simulated 
data.  This  was  done  also,  and  the  results  displayed  in 
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Figure  4.14.  The  root  aetn  square  error  by  the  nonptriaet- 
ric  method  is  a  factor  of  20  smaller  (1.021E-4  versus 
2.040E-3)  than  by  using  Shannon  maximum  entropy.  Thus,  the 

nonparaaetr ic  tool  cam  be  used  as  either  a  rigorous  numeri¬ 
cal  method  for  fiudimg  the  distributiom  of  a  fuuetiom  of 
multiple  random  variables,  or  as  a  direct  estimator  of 
distribution  functions  given  a  data  set. 

A  Comment  on  Sparse  Data  Sets 

The  problem  of  sparse  data  sets  for  the  inputs  still 
remains.  It  is  still  true,  that  if  one  is  attempting  to 
measure  the  failure  probability  for  a  highly  reliable  part, 
and  the  data  for  the  inputs  is  limited,  then  one  is  forced 
to  estimate  the  area  in  the  tail  of  some  distribution.  If 
one  is  using  parametric  models,  then  the  relative  error  can 
be  large  for  different  parametric  models.  For  example, 
Stancampiano  [32]  found  a  failure  probability  of  7E-10  +  5E- 
10  using  normal  fitting;  2E-31  +  2E-31  using  importance 
sampling,  and  2E-16  +  4E-16  using  Shannon  maximum  entropy. 
There  is  a  large  relative  variation  within  each  method  and 
between  the  methods,  even  though  the  absolute  value  is  close 
to  zero.  The  nonparame tr ic  estimate  for  the  same  problem  is 
zero . 

Tail  estimation  is  always  difficult.  Shapiro  and  Gross 
[12]  point  out  that  great  care  needs  to  be  exercised  when 
one  seeks  probability  values  that  are  outside  the  range  of 
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the  set  o£  dsts  obtained.  The  nse  of  an  erroneous  paramet¬ 
ric  model  can  be  costly  in  this  case.  The  advantage  of  the 
nonparametr ic  approach  taken  here  is  that  one  cannot  infer 
probabilities  far  beyond  the  range  of  the  data.  Protection 
against  unreasonable  inference  is  built  in. 

Further  research  in  applying  the  nonparametr ic  tool  to 
highly  reliable  assemblies  is  in  order.  However,  that  is 
not  a  problem  in  aircraft  nuclear  survivability.  As  men¬ 
tioned  earlier,  this  tail  issue  is  not  important  for  this 
application.  One  only  has  a  tail  problem  if  one  can 
guarantee  that  the  location  parameter  for  the  stress  vari¬ 
able  is  fixed  in  the  stress  space.  For  the  aircraft  nuclear 
survivability  problem  no  such  guarantee  can  be  given. 
Besides  this,  extremely  small  failure  probabilites  (sure- 
safe  regions)  and  extremely  large  failure  probabilites 
(sure-kill  regions)  are  not  the  regions  where  one  worries 
about  what  happens.  It  is  clear  what  happens.  It  is  the 
space  between  that  one  wonders  about,  and  any  techique  that 
can  accurately  find  the  shape  of  the  distributions  in  this 
region  is  a  valuable  tool.  Such  is  the  case  for  the  nonpar¬ 
ametr  ic  model. 

Summary 

To  summarize,  the  no np a r a m e t r  i  c  method  that  will  be 
used  to  find  the  survivability  of  airoraft  in  nuclear  weapon 
environments  has  been  described.  A  benchmark  problem  from 
the  reactor  risk  analysis  literature  has  been  worked,  and 
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the  remits  compared  to  more  standard  methods.  The  nonpara- 
metric  method  has  been  shown  to  be  numerically  efficient, 
and  nsefnl  as  a  stand-alone  numerical  algorithm  or  a  tool 
for  direct  application  on  existing  data  sets.  Its  use  on 
sparse  data  sets  has  been  mentioned.  Its  chief  advantage 
here  is  that  it  provides  protection  against  unreasonable 
extrapolation  beyond  the  information  contained  in  the  sparse 
data  set. 

In  the  next  two  Chapters  the  tool  just  described  and 
demonstrated  will  be  applied  to  the  problem  of  aircraft 
survivability  in  nuclear  weapon  environments. 


The  prescription  for  a  survivability  calculation  can 
now  be  spe c if ied . sue c inc t ly.  One  needs  to  first  construct  a 
reliability  model  of  the  system  as  a  whole.  For  each 
component  that  goes  into  the  model,  one  then  requires  a 
stress  and  a  strength  function  for  the  component.  For 
nuclear  survivability,  the  nuclear  weapon  effect  (in  this 
case  overpressure)  appears  in  the  stress  model.  If  the 
weapon  effect  is  a  random  variable,  then  the  applied  stress 
is  also.  The  strength  model  may  be  a  single  number,  that 
is,  a  sure-kill  (SK)  specification,  or  it  may  be  a 
functional  equation  as  in  Chapter  IV.  For  the  present 
Chapter,  a  single  number  will  suffice  provided  one  can  test 
enough  items  so  that  the  distribution  of  such  numbers  can  be 
determined.  This  distribution,  of  course,  is  the  strength 
distribution.  Having  obtained  both  distributions,  which  may 
vary  in  space  and  time,  the  failure  probability  is 
determine  d . 

Mission  Survivab i 1 i t v 

A  systems  approach  to  the  problem  of  nuclear  surviva¬ 
bility  analysis  must  begin  with  a  description  of  the 
intended  mission  and  the  threats  to  that  mission.  Every 


mission  may  be  thought  of  as  a  series  chain  of  mission 


phases.  Each  phase  must  be  executed  ixx  order  for  the  mis¬ 
sion  to  be  completed.  These  ideas  have  been  recognized 
[45].  An  illustration  of  the  base  escape  problem  (from  a 
paper  by  Bridgman  [45])  is  shown  in  Fignre  5.1. 

Daring  each  mission  phase,  the  weapon  system  may 
encounter  one  or  more  nuclear  bursts  at  particular  points  in 
space  and  time.  Thus  each  phase  may  be  modeled  as  a  chain 
of  burst  encounters.  Each  encounter  will  cause  various 
nuclear-induced  stresses  to  be  placed  on  the  weapon  system. 


weapon  induced  stress  some  subsystems  will  be  vulnerable  and 
some  will  not.  Those  subsystems  that  can  fail  must  be 
identified.  The  calculation  of  the  survivability  may  then 
proceed  as  suggested  by  Bridgman  [45].  His  flow  diagram 
from  Reference  [45]  is  shown  in  Figure  5.2. 

Two  of  the  most  difficult  parts  of  the  calculation  are 
shown  in  Figure  5.2  within  the  dashed  lines.  Regarding  the 
first  box,  'Find  Intensity  on  Target',  Bridgman  writes  "The 
calculation  of  free  field  intensities  ....  is  a  challeng¬ 
ing  business  involving  all  manner  of  calculation  on  the 
hydrodynamics  of  blast  waves  on  the  transport  of  neutrons, 
gamma  rays  and  x-rays  on  the  propagation  of  the  electromag¬ 
netic  pulse,  etc.  All  of  this  has  occupied  our  attention 
for  too  long  a  period.  Our  focus  here  is  to  go  beyond  the 
free  field  calculations  and  focus  on  the  next  block  of  the 
diagram — finding  the  probability  of  failure  .  .  ."  [46]. 


Figure  5.1.  Pictorial  of  Base  Escape 
(From  Reference  [45]) 


In  all  that  follows,  the  focus  of  this  dissertation 
will  he  limited  to  the  two  elements  of  the  block  diagram 
just  mentioned.  Actually,  one  can  generalize  those  two 
blocks  just  a  bit  more  by  considering  the  intensity  on  the 
target  as  a  random  variable  also.  The  scope  of  work  in  this 
Chapter  and  the  next  is  illustrated  in  Figure  5.3,  which 
expands  on  the  two  blocks  of  Figure  5.2.  Some  additional 
explanation  of  Figure  5.3  is  in  order. 

In  a  deterministic  model,  one  would  find  the  intensity 
of  the  stress  on  target  and  then  compare  that  stress  value 
to  the  strength  of  the  part.  If  the  intensity  exceeds  the 
strength  of  the  part,  then  the  part  fails.  Otherwise,  it 
does  not. 

If  one  uses  a  probabilistic  model,  then  one  finds,  not 
the  intensity  on  the  target,  but  the  statistical  distribu¬ 
tion  of  intensities  on  the  target  at  the  particular  range 
and  time  of  interest.  Also,  the  strength  of  the  item  is  not 
now  represented  by  a  single  number,  but  by  a  statistical 
distribution  of  numbers.  The  failure  probability,  which  is 
clearly  a  function  of  range  and  time  and  other  variables,  is 
now  given  by  the  interference  integral  of  Equation  3.3. 

The  central  problem  is  to  determine  the  distributions 
of  stress  and  strength  at  each  range  and  time  point  of 
interest.  There  are  three  ways  to  do  this,  and  each  is 
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First,  one  can  simply  claim  to  know  th<>  distributions, 
or  rather,  one  can  select  them  a  priori.  As  has  been  pre¬ 
viously  discussed,  selection  of  the  Dirac  delta  function  for 
the  stress  and  strength  PDF's  leads  directly  back  to  the 
cookie-cutter  way  of  doing  business.  One  might  select  the 
Dirac  delta  function  for  the  stress  PDF  and  a  lognormal  for 
the  strength  PDF.  This  is  precisely  what  is  illustrated  in 
Figure  5.2.  The  intensity  on  target  is  believed  known.  The 
part  fails  randomly.  Of  course,  the  disadvantage  of  this 
approach  is  that  it  may  be  difficult  to  defend  the  choice  of 
distributions  used. 

A  second  approach  is  the  straightforward  one  of  atmr- 
ing  the  quantities  desired.  The  stress  distribution  is 
found  by  measuring  the  stress  over  a  number  of  identical 
experiments  so  that  the  random  character  of  the  stress  is 
determined.  The  strength  distribution  is  determined  by 
testing  a  number  of  items  to  complete  failure.  This  method 
may  be  termed  a  macrodata  approach,  since  the  information 
obtained  is  directly  applied  to  the  failure  calculation  of 
the  system  as  a  whole.  The  disadvantage  is  that  the  testing 
requirements  may  be  impossible. 

Finally,  one  may  infer  the  quantities  desired.  A 
deterministic  model  of  the  stress  is  formulated  as  in,  for 
example,  the  bending  moment  on  a  wing  of  an  aircraft.  This 
model  may  depend  on  other  variables  whose  distributions  are 
known.  The  statistical  distribution  of  the  stress  (bending 


moment)  is  then  found  by  propagating  the  statistical  varia¬ 
tion  of  all  the  input  variables  into  the  output  variable. 
The  strength  distribution  is  found  similarly.  A  good  exam¬ 
ple  of  an  inferred  strength  distribution  is  given  by  Equa¬ 
tion  4.11  which  predicts  the  burst  pressure  for  a  reactor 
pressure  vessel  [32].  One  might  call  this  a  microdata 
approach,  since  the  information  collected  is  not  applied 
directly  to  the  system  of  interest  but  must  be  propagated 
into  a  larger  model  before  one  can  calculate  the  failure 
probability.  Although  this  approach  is  more  rigorous  and 
useful  than  the  others  considered,  it  has  one  potentially 
serious  disadvantage.  It  depends  critically  on  an  accurate 
model  of  the  processes  being  considered.  Unlike  direct 
testing,  one  is  calculating  the  probability  of  failure  given 
the  modeling  is  accurate.  Catastrophic  events  may  well 
occur  if  the  modeling  is  incomplete  or  inaccurate. 

The  central  thrust  of  this  Chapter  will  be  to  determine 
the  stress  and  strength  distributions  for  selected  aircraft 
components  vulnerable  to  the  effects  of  overpressure  and 
blast  from  a  nuclear  weapon.  All  of  the  above  approaches 
will  be  considered,  including  mixtures  of  the  approaches. 

Stress  Distr ibut ions  in  a  1  Kiloton  Blast  Environment 

The  overpressure  environment  from  a  nuclear  burst  in  a 
homogeneous  atmosphere  may  be  described  in  three  separate 
regions.  These  are  the  free  air  region,  the  transition 
region,  and  the  mach  stem  region.  In  the  free  air  region. 


no  ground  effects  perturb  the  blast  live.  In  tbe  each  stem 
region,  the  reflected  blast  wave  has  caught  up  and  construc¬ 
tively  interfered  with  the  primary  wave  so  that  a  region  of 
intense  overpressure  exists  near  the  ground.  The  transition 
region  is  the  region  between  the  free  air  and  mach  stem 
regimes,  typically  in  the  vicinity  above  the  triple  point 
[45]  . 

The  overpressure  environment  in  the  free  air  region  is 
usually  described  by  the  Air  Force  Weapons  Laboratory  (AFWL) 
1  Kiloton  (KT)  standard  [74],  or  the  more  recent  DNA  stan¬ 
dard  [73].  Research  has  been  conducted  to  further  refine 
this  standard  [47,51].  Reference  [47]  includes  the  95% 
confidence  intervals  about  the  nominal  pressure  versus  range 
curve.  This  is  shown  in  Figure  5.4,  which  is  a  reproduction 
from  that  Reference.  This  information  can  be  analyzed  to 
construct  the  statistical  distribution  of  overpressure  as  a 
function  of  range  from  the  nuclear  burst. 

Polynomial  fits  exist  for  the  nominal  pressure  versus 
range  curve  and  are  available  in  the  literature  [ 66 ]. 
Unfortunately,  Reference  [47]  does  not  give  the  polynomial 
fits  to  either  of  the  bounding  curves  that  represent  the  95% 
confidence  interval.  Since  this  information  is  lacking,  a 
direct  analysis  of  the  graph  of  Figure  5.4  was  made. 
Although  the  dispersion  was  found  to  increase  slightly  with 
range,  an  average  value  was  used  in  the  following  analysis. 


Overpressure  P(psig) 


Figure  5.4  is  a  log-log  plot  o£  p,  the  overpressure, 
versus  r,  the  range  from  the  burst.  If  it  is  assumed  that  a 
normal  regression  of  In  p  on  In  r  has  been  done,  then,  for  a 
fixed  value  of  r,  the  random  variable  P  is  log-normal  with 
location  parameter  and  scale  parameter  Pp.  Note  that  the 
lognormal  distribution  follows  naturally  from  the  standard 
assumption  that  the  distribution  of  the  residuals  in  the  log 
space  is  normal.  The  location  parameter  is  given  by  the 
regression  equation,  while  the  scale  parameter  is  assumed 
constant.  For  any  range  r  then,  measured  in  meters,  the 
distribution  of  P  (in  Pascals)  is  lognormal  with  parameters 
given  by: 

a  (r)=.19(ln(r/1000))2-l.Sln(r/1000)+8.738  (5.1) 

P 

pp  (  r )  =.l 7  (5.2) 

The  stress  distribution  far  particular  components  will 
take  its  random  nature  from  the  random  character  of  the 
basic  overpressure  variable.  In  one  case,  this  results  in  a 
simple  analytic  problem  for  the  failure  probability,  but  in 
most  cases  one  must  find  the  distribution  of  functions  of 
the  overpressure  variable. 

Strength  Distributions  for  A_i r.c.r j»X_t  Components 

With  the  fundamental  random  variable  input  for  the 
stress  function  identified,  one  can  now  consider  strength 
functions  and  strength  distributions.  Three  approaches  will 
be  examined,  all  of  which  have  appeared  in  the  literature  at 
one  time  or  another.  These  include  cookie-cutter  techniques 


[54],  probabilistic  methods  with  an  a  priori  lognormal 
assumption  [45],  and  actual  use  of  long-term  test  data 
[22.20]  . 

Cookie-Cutter  ]) Al.£ r Aiu£.i o.n .  For  all  of  the  work  in 
this  Chapter,  the  strength  model  is  given  by  the  simple 
equation : 

S-SK  (5.3) 

That  is,  the  single  sure-kill  value  SK  characterizes  the 
strength  of  the  component.  For  the  c o ok i e- cu t t e r  case,  as 
noted  in  Chapter  III,  the  strength  PDF  is  then  given  by: 

fs(s)~3(s-SK)  (5.4) 

and  the  strength  CDF  by: 

Fs(s)=0  v  s  < SK  (5.5) 

Fs(s)=*l  v  s2.SE  (5.6) 

The  parameter  for  the  Dirac  delta  distribution  is  the  sure- 
kill  specification,  given  by  sources  such,  as  Reference  [5^]. 
Some  references  [42]  use  an  unbiased  cookie-cutter  speci¬ 
fication  which  uses  the  average  of  the  sure-safe  (SS)  and  SK 
specifications  as  the  parameter  of  the  Dirac  delta  function. 

These  techniques,  which  remain  in  practice  [56],  result 
in  cookie-cutter  damage  distributions  with  range,  provided 
the  stress  distribution  is  also  a  Dirac  delta  function. 

A  Pr.io.rjL  Lo^noraal  Distribution.  A  second  common 
approach  is  the  assumption  of  a  lognormal  strength 
distribution,  with  a  concurrent  estimation  of  its  parameters 
by  some  method  [45,42,50,71]. 


Bridgman,  for  example  [45],  assigns  the  SS  and  SE 
spe  o  i  f  i  c  a  t  i  on  s  from  Reference  [54]  as  the  2nd  and  98th 
percentiles,  respectively,  of  the  strength  distribution  of  a 
specified  component.  Denoting  SS  and  SE  by  S  Q2  and  S  g 
respectively,  the  parameters  of  the  lognormal  are  found  from 
the  equations: 

“S={la<S.02S.98>}/2  (5.7) 

f  ln(  s.9g/ S.o2^  *  ^2z.9  8*  (5.8) 

where  z  is  the  98th  percentile  of  the  standard  normal 
variate  and  is  approximately  2.054. 

This  approach  satisfies  qualitatively  the  concepts  of 
SS  and  SE  specifications,  while  at  the  same  time  providing  a 
continuous  probability  of  failure  with  range  from  the 
nuclear  burst. 

Strength  Distributions  From  Test  Data.  The  cookie- 
cutter  approach  does  not  take  into  account  the  probabilistic 
nature  of  failure  phenomenon,  while  the  alternate  approaches 
of  Bridgman  [45]  and  others  [42],  though  probabilistic,  lack 
foundation  in  the  actual  database,  and  do  not  treat  the 
uncertain  nature  of  the  applied  stress  environment.  Conse¬ 
quently,  the  Department  of  De.fense  (DOD)  database  for 
nuclear  gust  and  overpressure  effects  on  aircraft  was  care¬ 
fully  surveyed  to  find  aircraft  component  failure  data.  The 
results  of  that  data  search  are  presented  at  this  time. 


POD  Database  Fog  Blast  Effects.  The  DOD  database 
for  blast  effects  on  aircraft  and  aircraft  structures  falls 
into  two  broad  classes — simulation  testing  and  nuclear  test¬ 
ing.  The  Defense  Atomic  Support  Information  and  Analysis 
Center  (DASIAC)  has  compiled  complete  bibliographies  on  both 
simulation  and  nuclear  testing. 

Simulation  testing  is  largely  carried  out  by  the 
Defense  Nuclear  Agency  (DNA),  though  others  are  involved  in 
it  also.  A  complete  list  of  DNA  and  non-DNA  reports  on 
simulation  testing  may  be  found  in  Reference  [35].  The 
purpose  of  simulation  testing  is  to  compare  theoretical 
response  predictions  with  actual  response  observations. 
There  is  little  failure  information  in  this  database,  since 
the  objective  of  simulation  testing  is  not  to  test  to 
destruction.  Failures  do  occur  sometimes  in  simulation 
experiments,  but  these  are  almost  always  unintentional.  An 
example  of  this  can  be  found  in  Reference  [60]  where  the 
tail  section  of  an  A-4C  was  destroyed  in  a  simulation  of  a  1 
KT  yield  device. 

Nuclear  testing  has  been  carried  out  both  above  ground, 
prior  to  1963,  and  by  underground  tests.  Reference  [35] 
also  lists  reports  compiled  from  nuclear  tests,  most  of  them 
above  ground.  In  addition,  DASIAC  has  published  two  unclas¬ 
sified  summaries  of  operations  REDWING  [70]  and  PLDMBBOB 
[69].  Besides  these,  DASIAC  has  also  published  a  complete 
compendium  of  blast  and  thermal  effects  on  aircraft  gleaned 


from  a  number  of  above  ground  nuclear  teats  [64].  Table  V.l 
above  tbe  operations  and  shots  that  were  examined.  Some  of 
these  tests  involved  in-flight  aircraft.  Since  these  were 
manned,  great  care  was  taken  to  be  sure  that  applied  stress 
levels  were  well  below  limit  loads  of  the  components.  Typi¬ 
cal  response  data  for  these  is  given  in  Table  V.2  taken  from 
Reference  [70].  The  response  data  is  displayed  as  a  frac¬ 
tion  of  design  limit  load. 

Although  no  failure  information  was  found  for  in-flight 
aircraft  response  to  nuclear  overpressure  effects,  much  is 
available  for  parked  aircraft.  Unfortunately,  this  data 
consists  mostly  of  verbal  descriptions  of  damage  to  a  large 
number  of  different  components.  In  many  cases  the  exact 
failure  mechanism  is  not  known  (e.g.--'The  canopy  was 
cracked  and  badly  burned.')  An  attempt  was  made  to  quantify 
the  information  by  constructing  man-hour s-to-repa ir  versus 
overpressure  plots,  but  these  are  difficult  to  extrapolate 
and  remain  classified.  Consequently,  the  nuclear  tests 
database  does  not  contain  the  information  needed  to  deter¬ 
mine  statistical  failure  distributions  for  aircraft  compo¬ 
nents. 

Even  though  the  DOD  database  did  not  yield  information 
that  would  allow  reliability  theory  to  be  used,  useful 
information  was  obtained.  If  drone  aircraft  in-flight  tests 
had  been  conducted  at  higher  levels,  one  might  have  observed 
failures  or  a  number  of  failures.  If  the  number  of  failures 
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TABLE  V.l 


NUCLEAR  EVENTS  INVOLVING  VULNERABILITY  TESTING 


Vainer ab il ity 

Operation 

Nuclear  Shot 

A/C  In-Flight 

UPSHOT-KNOTHOLE 

HARRY 

NANCY 

SIMON 

No  Failures 

Observed — Max  Load 
of  1.1  *  Limit  Load 

REDWING 

DAKOTA 

NAVAJO 

HURON 

APACHE 

HARDTACK 

KOA 

WALNUT 

DOMINIC 

AZTEC 

ENCINO 

YESO 

A/C  Structures+ 

UPSHOT-KNOTHOLE 

ENCORE 

GRABLE 

+  Yield  Points 
Exceeded  in  2  Cases 


?.1« 


TABLE  V .  2 


REPRESENTATIVE  RESPONSE  DATA  (From  Reference  [70]) 


T0  Is  time 

SUfWARY  DATA  ON 

of  detonation. 

PROJECT  5.3 
Tsa  Is  time 

of  shock  arrival . 

Abso 1 ute  Hor 1 zon  ta  1  Range 

Horizontal  Range 

Radiant 

Peak 

Percent  of 

Shot 

Altitude 

at  T0 

Tsa 

Exposure 

Overpressure 

Design  Limit 

feet 

feet 

feet 

ca 1 /cm2 

psl 

Cherokee 

34,000 

47,785 

139,571 

52  Wing 

ZunI 

14,000 

27,000 

97,760 

34.0 

0.400 

47  Wing 

Flathead 

16,000 

17,800 

59,100 

45  Wing 

Dakota 

16,000 

13,100 

35,050 

88  Wing 

Apache 

8,000 

23,500 

60,500 

48  Wing 

Navajo 

— 

— 

— 

46.2  Wing 

Huron 

9,894 

8,768 

23,386 

110  Wing 

Tewa 

19,000 

27,250 

65,750 

0.900 

65  Wing 

T0  Is  time 

SUMMARY  DATA  ON 

of  detonation. 

PROJECT  5.4 
Tsa  Is  time 

of  shock  arrival. 

Absolute  Horizontal  Range 

Horizontal  Range 

Radiant 

Peak 

Percent  of 

Shot 

Altitude 

at  TQ 

«t  Tsa 

Exposure 

Overpressure 

Design  Limit 

feet 

feet 

feet 

cal /cm2 

psl 

Lacrosse 

13,700 

6,750 

28,200 

1.17 

0.283 

35  Wing 

Zun  1 

16,900 

34,000 

— 

10.55 

— 

Erie 

10,450 

3,829 

14,000 

65  Wing 

Flathead 

25,700 

13,466 

44,659 

50  Wing 

Inca 

9,815 

2,624 

11,572 

— 

Dakota 

17,650 

16,690 

43,630 

56  Wing 

Apache 

10,200 

28,516 

45,375 

35  Wing 

Huron 

16,200 

10,000 

31,493 

40  Wing 

TQ  Is  time 

SUMMARY  DATA  ON 

of  detonation. 

PROJECT  5.5  CAPABILITIES  AIRCRAFT  (F-84F) 

Tsa  Is  time  of  shock  arrival. 

Absolute  Horizontal  Range 

Horizontal  Range 

Radiant 

Peak 

Percent  of 

Shot 

Altitude 

at  T0 

at  Tsa 

Exposure 

Overpressure 

Design  Limit 

feet 

feet 

feet 

cal /cm2 

psl 

Lacrosse 

15,200 

8,640 

600 

1.65 

0.80 

51  Wing 

Flathead 

21,000 

900 

14,935 

40  Wing 

Dakota 

20,374 

300 

11,065 

115  Wing 

Mohawk 

19,920 

3,729 

19,700 

65  Wing 

Apache 

31,614 

18,614 

64,300 

— 

Navajo 

16,865 

25,200 

57,400 

10  Wing 

had  bean  cataloged  as  a  function  of  the  ratio  of  stress  to 


design  allowable  stress,  a  generic  strength  distribution 
might  be  constructed  which  should  apply  to  any  component 
that  is  produced  by  the  same  design  process.  Such  data  was 
indeed  found  in  non-DOD  publications,  as  discussed  below. 

St.  .a  t.  i.  c.  Te.s„ts,  On  A  i.r.c  ra„f  t,.  A  computer  search  of 
the  DIALOG  (a  commercial  literature  search  service)  data¬ 
bases  located  two  separate  compilations  of  static  test 
results  on  aircraft.  The  first  is  a  compilation  by  Jablecki 
[26]  done  in  1955,  and  based  on  static  test  results  con¬ 
ducted  at  V r i g h t - P a t t e r s o n  AFB  (VPAFB)  from  1940-49. 
Jablecki  did  not  actually  do  a  statistical  analysis.  This 
was  later  done  both  by  Chenoweth  [22]  and  by  Freudenthal  and 
Wang  [24].  A  second  set  of  static  test  data  was  obtained 
and  analyzed  by  Chenoweth  in  1972  [20].  This  was  a  compila¬ 
tion  of  static  tests  conducted  at  WPAFB  from  1950-1970.  In 
both  cases,  the  data  involved  static  loadings  of  the  wings, 
fuselage,  horizontal  stabilizer,  vertical  stabilizer,  and 
landing  gear.  The  data  analysis  was  done  in  a  normalized 
fashion  by  constructing  the  failure  probability  as  a  func¬ 
tion  of  the  ratio  of  applied  stress  to  design  ultimate 
stress.  These  analyses  allow  one  to  construct  the  failure 
distribution  for  an  aircraft  component  as  a  function  of  the 
design  parameters.  Both  of  these  data  sets  have  been  used 
in  estimating  the  reliability  of  aircraft  structures  to 
thunderstorm  gust  loadings  [24,21]  and  aircraft  maneuver 


loadings  [21].  A  more  detailed  discussion  of  these  findings 
follows . 

Large  Aircraft  Structure  s.  Chenoweth  fit  the 
Jablecki  data  for  the  larger  structures,  the  wings  and 
fuselage,  to  a  lognormal  failure  distribution  with  the  ran¬ 
dom  variable  being 

I  =  S/SUL  (5.9) 

where  the  random  variable  S  is  the  stress  on  the  component 
that  causes  it  to  fail  and  Sp^  is  tlie  ultimate  load  design 
stress  for  the  component.  Consequently,  the  component 

failure  distribution  is  given  by 

Fx(x)=d>((lnx-ax)/px))  ▼  x>0  (5.10) 

where  <t>  ( z )  is  the  standard  normal  cumulative  distribution 

function  (CDF)  evaluated  at  z,  and  given  by 

<p(z)=Jl'D(2n)~1/2exp(-t2 /2)dt  (5.11) 

The  variable  x  is  an  arbitrary  value  of  the  random  variable 
X,  <*x  is  the  location  parameter  for  the  lognormal  density 
function,  and  0X  is  the  scale  parameter  for  the  lognormal 
density  function.  If  one  desires  the  distribution  of  the 
random  variable  S  it  can  be  shown  that  given  X  is  lognormal 


with  parameters  ax  and  fJ  x  ,  then  S  is  lognormal  with 


parameters  <*s  and  Pg  given  by 
aS~  °x  +  1  nS0L 

PS  =  Px 


(5.12) 

(5.13) 


Small  Aircraft  Structure  s.  Chenoweth  fit  the 
data  from  the  smaller  structures  (vertical  stabilizer. 


horizontal  stabilizer,  and  landing  gear)  to  an  asymptotic 
Weibull  distribution  function.  This  is  just  a  power 
function,  and  the  distribution  of  X  in  this  case  is  given  by 
Fx(x)=0  v  xi.0  (5.14) 

FI(x)  =  (x/ax)1^bx  v  0£xiax  (5.15) 

Fx(x)=l  v  x2ax  (5.16) 

In  the  above  equations  ax  and  bx  are  the  parameters  for 
the  power  function  distribution.  If  the  distribution  of  S 
is  desired  it  can  be  shown  that,  given  X  has  a  power 
function  distribution  with  parameters  ax  and  bx  then  S  has  a 
power  function  distribution  with  parameters  a§  and  bg  given 
by 

aS“axSUL  *  (5.17) 

bs=bx  (5.18) 

The  1972  data  analysis  is  not  broken  out  into  large  and 
small  structures  as  was  done  for  the  Jablecki  data.  The 

data  is  instead  lumped  together  and  fitted  to  the  power 
function  model.  Chenoweth  shows  [20]  both  sets  of  data 
fitted  to  this  model,  and  notes  that  a  slight  improvement 
can  be  seen  in  the  reliability  of  the  design  process  from 
1945  to  1960.  The  improvement  is  shown  in  Figure  5.5, 
taken  from  Reference  [20].  For  the  more  recent  data  set, 
the  mean  has  improved  by  about  10%  and  the  coefficient  of 
variation  (ratio  of  standard  deviation  to  mean)  has  been 
reduced  by  about  33%. 
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Cite  no  we  th's  date  is  somewhat  aurprising.  It  aeeaa  to 
iaply  that  at  noainal  liait  loads  (2/3  of  nltiaate)  the 
failure  probability  is  significantly  high — on  the  order  of 
10%.  Chenoweth  apparently  insists  on  the  validity  of  the 
data.  Quoting  froa  Reference  [21],  he  aays,  "it  ia  ahown 
that  a  'no  static  test1  policy  on  aajor  aircraft  strnctnral 
coaponents  yields  extreaely  high  probability  of  failnre  or 
unreliability."  No  foraal  rebuttal  to  this  was  found  in  the 
literature.  Chenoweth  goes  on  in  Reference  [21]  to  further 
show,  based  on  the  data,  the  vulnerability  of  aircraft  to 
thunderstorm  gust  loads  and  aanenver  loads.  In  both  cases, 
the  resulting  failure  probabilities  range  froa  IE-4  to  IE-2. 
The  following  additional  observations  should  also  be  noted. 

First,  the  data  as  shown  in  Figure  5.5  is,  at  sost, 
half  the  story.  The  distribution  of  stresses  or  loads  can 
be  as  important,  or  perhaps  more  important,  than  the 
strength  distribution  in  determining  the  actual  probability 
of  failure.  For  example,  if  the  ultimate  load  design  had 
been  chosen  in  such  a  way  that  the  limit  load  value  (2/3  of 
ultimate,  say)  was  extremely  improbable,  then  the  aircraft 
could  be  perfectly  hard.  It  is  the  interference  of  the 
stress  distribution  with  the  strength  distribution  that 
governs  the  failure  probability,  and  not  either  acting 
a  1  one . 

As  a  second  consideration,  Chenoweth  says  [21]  the  data 


is  first  time  static  test  data.  If  these  tests  resulted  in 


redesign  of  the  components,  then  the  resulting  product  would 
of  course  he  herder  than  the  first  time  tests  indicate. 

Third,  the  data  is  admittedly  old.  The  most  recent 
data  set  is  a  mixture  of  aircraft  from  1950  to  1970.  The 
methods  of  aircraft  design  and  hardness  verification  and 
validation  have  undoubtedly  changed  since  the  data  was 
obtained,  so  that  one  should  not  infer  that  it  applies  to 
modern  combat  airoraft.  Because  of  this,  the  resulting 
calculations  in  this  Chapter  should  be  considered  worst 
case  results. 

Having  said  all  this,  the  data  is  still  a  valuable 
collection  of  information.  Here  at  least  is  some 
justification  for  choosing  parametric  failure  distributions. 
The  following  calculations  illustrate  the  use  of  such 
information. 

Component  Fa ilure  Probab i 1 i t ie  s 

Flilliii  The  general  failure 

probability  is,  as  stated  earlier,  given  by  the  reliability 
interference  integral.  That  is, 

”f  s  (s)Fs(*)d«  (5.19) 

Model  1.  The  result  depends  on  the  time  dependent 
properties  of  both  the  strength  and  stress  distributions. 
The  stress  function  is  presumed  given  by 

s«P  (5.20) 


where  P  is  just  the  random  overpressure  variable, 
lognormally  distributed  with  parameters  given  by  Equations 


5.1  and  5.2.  For  this  case,  daaoted  as  Model  1.  the  cookie- 
cutter  strength  distribntion  is  assumed,  reducing  Equation 
5.19  to  the  re  suit: 

pf«l-F#(SK)  (5.21) 
where  SK  is  1.085E5  Pa  (15.75  psi).  Since  the  stress 
distribution  is  a  lognormal,  the  results  for  this  case  have 
an  analytical  form  as  a  function  of  range.  Explicitly, 

pf(r)-l-+((lnSK-ap(r))/pp)  (5.22) 
where  <t>(z)  is  given  by  Equation  5.11. 

The  results  of  Model  1  are  shown  in  both  the  stress 
space  (Figure  5.6)  and  the  range  space  (Figure  5.7).  One 
notes  that  the  use  of  a  probabilistic  model  for  the  stress 
is  sufficient  to  provide  a  continuous  probability  of  failure 
in  the  range  space.  The  tails  in  Figure  5.7  for  Model  1  are 
due  to  the  variation  in  the  applied  overpressure  only,  since 
the  strength  distribution  is  a  step  function.  The  median 
failure  range  for  Model  1  corresponds  to  the  SK  stress 
spe  c if ica t ion . 

Mod9l  2.  Keeping  the  stress  distribution  the 
same,  if  the  lognormal  assumption  is  made  for  the  strength 
distribution,  with  parameters  calculated  from  Equations  5.7 
and  5.8.  and  SS=6891  Pa  (1  psi),  and  SK-1.085E5  Pa  (15.75 
psi),  then  the  failure  probability  may  be  written  as: 

pf-Pr{SiP}-y]®fp(p)Fs(p)dp  (5.23) 
where  S  is  the  strength  random  variable,  P  the  stress  random 
variable  (the  applied  overpressure),  f»(p)  the  probability 


density  function  of  P,  and  Fg(p)  the  CDF  of  S  evaluated  at 
p.  Equation  S.23  could  also  be  written  as 


pf«Pr{S/PilJ  (5.24) 

or 

p£-Pr(ln(S/P)iO)  (5.25) 

In  the  case  of  Equation  5.25.  letting  the  random 
variable  Y,  sometimes  oalled  the  log  safety  factor.  be 
defined  by 

Y»ln(S/P)-lnS-lnP  (5.26) 

it  is  seen  that  Y  is  consequently  normally  distributed  with 
mean  Py  and  standard  deviation  g given  by 

•*y"ttS“ap  (5.27) 

dy-(Ps2  +  Pp2)1/2  (5.28) 

Consequently,  the  failure  probability  for  the  fuselage 
can  be  shown  to  be 

pf-4>(-|i  y/tfy)  (5.29) 

Since  the  parameter  <*p  is  a  function  of  distance  from 
the  burst,  p£  as  given  by  Equation  5.29  is  also  a  function 
of  distance  from  the  burst.  The  results  for  Model  2  are 
also  shown  in  Figures  5.6  and  5.7. 

M 9 j p !•  Finally,  the  actual  data  for  fuselage 
failures  may  be  used  to  form  the  strength  distribution  while 
retaining  the  previous  model  for  the  stress  distribution. 
Since  the  large  structures  were  found  to  have  lognormal 
strength  distributions,  the  general  form  for  the  failure 
probability  is  given  again  by  Equations  5.27  through  5.29. 


However,  in  this  ease  the  piraaeteri  <*g  end  pg  ere 
calculated  froa  Equations  5.12  end  5.13  where  the  pereaeter 
S|jl,  i*  related  to  the  anre-aefe  specif ieetion  by: 


Sul-1.5(SS)  (5.30) 

end  the  peraaeters  SS,  <*x,  end  0Z  ere  |iren  by: 

SS-6891  Pe  (1  psi)  (5.31) 

ax«.1481  (5.32) 

Px«.3159  (5.33) 


The  vulnerability  of  the  fuselage  for  this  case  is  also 
illustrated  in  Figures  5.6  and  5.7.  Model  3  (based  on  the 
Chenoweth  data)  aeeas  to  iaply  a  auch  softer  fuselage  than 
either  Model  1  or  Model  2  would  suggest.  The  failure  proba¬ 
bility  rises  steeply  with  increasing  stress,  with  a  aedian 
failure  level  at  about  11.7B3  Pa  (1.7  psi).  Nevertheless, 
all  three  distributions  agree  qualitatively  in  the  SS  and  SK 
regiaes . 

Not  surprisingly,  the  softness  of  the  fuselage  as  a 
function  of  range  froa  the  1  KT  standard  is  clearly 
deaonstrated  in  Figure  5.7.  If  Model  3,  based  on  the 
Chenoweth  data,  is  accepted,  the  sure-safe  range  (denoted  by 
RS — the  distance  where  the  failure  probability  is  less  than 
or  equal  to  .02)  is  at  1160  aeters,  where  the  aedian 
overpressure  is  expected  to  be  5012  pascals  (.73  psi).  This 
stress  value  is  about  27%  smaller  than  that  postulated  by 
2048,  but  is  within  the  error  bounds  given  by  2048  [54]. 


Model 


This  soft  fuselage  result  qualitatively  agrees  vith 
that  of  the  work  documented  in  Reference  [60].  There  it  vas 
noted  that  the  2048  Method  1  technique  for  overpressure  vas 
severely  unconservative  .  In  a  test  of  the  A-4C,  the  SK 
condition  vas  reached  at  a  farther  range  from  a  1  IT  simula¬ 
tion  than  vas  expected.  In  fact*  the  catastrophic  failure 
occurred  at  311  meters  from  ground  zero.  From  Figure  5.7 
(keeping  in  mind  that  our  SS  and  SK  specifications  are 
somevhat  different  from  that  of  Reference  [60])  one  sees 
that  the  cookie-cutter  prediction  is  indeed  one  of  no  fail¬ 
ures,  vhile  Bridgman's  model  (Model  2)  vould  predict  a 
failure  probability  of  .76  at  that  point*  and  the  actual 
test  data  vould  predict  the  fuselage  to  be  sure-killed. 

Ag  ain,  the  caveats  previously  mentioned  should  be 
remembered.  In  addition  to  those,  there  can  still  be  an 
ambiguity  regarding  the  uses  of  the  terms  SS  and  SK. 
Although  Reference  [68]  clearly  equates  a  SS  condition  to  a 
limit  load  condition,  it  is  not  clear  vhat  the  factor  of 
safety  is  betveen  SK  and  SS  specifications.  Indeed,  for  the 
fuselage  specs  it  may  be  15.75  —  the  ratio  of  SK  to  SS.  If 
that  is  the  case,  then  the  ultimate  load  should  perhaps  be 
given  directly  by  15.75  psi  and  not  1.5  times  the  SS 
specification.  This  vould  have  a  dramatic  effect  on  the 
strength  distribution,  moving  the  median  failure  point  to  a 
higher  level. 


... . . . « 
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Wing 


Mode.  For  the  case  of 


the  fuselage,  the  presumption  that  the  fuselage  responds 
directly  to  overpressure  led  to  simple  computational 
results.  If  instead  a  more  complex  response  (i.e,  stress) 
model  had  been  chosen,  no  such  simple  results  would  have 
resulted.  This  is  the  case  for  even  elementary  models  of 
wing-loading . 

The  general  expression  for  the  failure  probability  as 
shown  in  Equation  5.19  is  still  valid,  but  now  has  no 
closed-form  solution.  Although  the  distribution  of  P 
remains  the  same  as  in  the  fuselage  analysis,  the  stress 
variable  a  is  no  longer  the  direct  overpressure  variable. 
Ije  stress  (using  a  simple  fundamental  mode  analysis  as  in 
Reference  [45])  is  now  the  load  factor  (in  g's — multiples  of 


weight)  given  by: 

nL=.5pv2CLA/W  (5.34) 

where  p  is  the  air  density  behind  the  shock  front  (kg/m2),  v 
is  the  lifting  wind  velocity  (m/sec),  CL  is  the  coefficient 
of  lift  (dimensionless),  A  is  the  area  of  the  wing  (square 
meters)  and  V  is  the  aircraft  weight  (Nt).  At  every  range 
point  r  the  distribution  of  the  load  factor  is  now  the 
required  stress  distribution.  The  problem  is  analytically 
intractable,  since  p,  v,  and  even  in  Equation  5.34  depend 
nonlinearly  on  the  overpressure. 

As  previously  discussed,  the  general  problem  of 
propagating  a  random  variable  through  a  complicated  response 


tits  staple  strength  or  response  function  models  justify. 
Fortunately,  nonpar an e t r i c  estimation  provides  a  powerful 
tool  [14, Appendix  A, Appendix  C]. 


The  distribution  of  load  factor  can  be  constructed  at 
any  range  r  by  using  the  technique  illustrated  in  Chapter 
IV.  This  is  done  by  drawing  a  stylized  saaple  froa  the 
assuaed  known  distribution  of  P  and  propagating  these  saaple 
values  through  the  response  equation  to  generate  a  load 
factor  sample.  This  saaple  can  then  be  used  to  soa'tuct 
the  distribution  [14, Appendix  A, Appendix  C]. 

With  this  tool  in  hand,  the  distribution  of  a  (load 
factor)  can  be  deterained  without  a  priori  distributional 
assumptions,  and  without  large-saaple  Monte  Carlo  analysis. 
Equation  5.19  can  then  be  numerically  integrated  at  each 
point  r  where  the  failure  probability  is  desired. 

For  the  case  of  wing-loading,  the  stress  model  is 
described  by  the  following,  where  bold-faoed  variables 
indicate  randoa  variables.  For  every  value  of  P  drawn  at 
range  r, 

(1)  Determine  the  magnitude  of  the  peak  wind  velocity, 
n^,  behind  the  shock  front  (from  Rankine-Hugoniot). 

a  / 


I 
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(A  complete  description  of  these  variables  may  be  found 
in  the  section  on  noaenclature.) 

(2)  Find  the  components  of  in  the  aircraft  frame  of 
reference . 

»l.x“(»l/*r>yb  (5.36) 

nj  ay=(»1/sr)(xt-xb)  (5.37) 

,lsia(,l/,r)di  (5.38) 

(3)  Determine  the  initial  attack  angle  (straight  and 
level  flight  assumed). 


clo“2W/(Povo  A) 


a0a‘(CL0“*356)/1*243 


rx  lax 
*ry— lay-v0 


(5.39) 

(5.40) 

(4)  Determine  the  resultant  vind  vector  in  the 

aircraft  frame  of  reference. 

(5.41) 

(5.42) 

*ri=*la*  (5*43) 

(5)  Resolve  the  resultant  vind  vector  into  its  lifting 

and  non-lifting  components. 

(5.44) 

(5.45) 

(6)  Determine  the  angle  between  the  lift  vind  and  the 
aircraft  chord  vector. 

«oadel=(cos(«o)llry+*in(aro^®rx  (5.46) 

delsAcos(eosdel)  (5.47) 

«=+ (n-d  e  1 )  (5.48) 


*pK*rx 

Y=(U  2  +  U  2)1/2 

v  ry  ii  ' 


(7)  Determine  the  coefficient  of  lift 


CL“1 .243<*+.3  56  (5.49) 

(8)  Determine  the  sir  density  behind  the  shock  front. 


P=P0{7+6P/pa}/(7+P/pa}  (5.50) 

(9)  Determine  the  load  factor  (a  sample  value  of  a) 
from  Equation  5.34,  i.e, 

nL=.5pY2CLA/W  (5.51) 

Mode  1  1.  With  the  stress  distribution  determined 
nonparametrically,  the  strength  distribution  is  now 
required.  Again  one  may  use  the  simple  equation, 

S  =  nsk  (5.52) 

For  Model  1  a  c ookie- cu t t e r  strength  distribution  is  assumed 
so  that  the  strength  variable  has  the  PDF, 

f  s(  s)  =  6(  s-nsk)  (5.53) 

where  nJk  is  taken  from  Reference  [54]  as  6.75  g's  for 
loading  from  below. 

As  in  the  fuselage  case,  the  results  for  the  failure 
probability  reduce  to 

pf  =  1-Fg(nsk)  (5.54) 

where  the  CDF  of  the  stress  distribution  is  now  found 
nonparametrically. 

The  results  for  loading  from  below  are  displayed  in 
Figures  5.8  and  5.9.  The  failure  probability  as  a  function 
of  applied  stress  is  shown  in  Figure  5.8,  while  the 
corresponding  failure  probability  with  range  is  shown  in 
Figure  5.9.  The  results  for  Model  1  are  shown  by  the  solid 


1  ine  s 


Model 


Mo4e 1  2 


The  stress  distribution  remains  the 


seme,  but  the  lognormal  assumption  is  now  made  for  the 
strength  distribution,  with  the  parameters  calculated  from 
Equations  5.7  and  5.8,  and  with  SS  and  SK  given  by  3  and 
6.75  g's  respectively  (Method  1  of  Reference  [54]). 

Since  the  stress  distribution  is  known  nonpar ame tr ic- 
ally,  the  integral  in  Equation  5.19  must  be  solved  numeric¬ 
ally.  The  method  for  doing  this  is  described  in  Appendix  D. 
The  results  using  the  lognormal  assumption  are  shown  in 
Figures  5.8  and  5.9  also. 

£.  Finally,  the  strength  distribution 

derived  from  the  test  data  is  given  as  a  lognormal,  with 
parameters  and  Pg  calculated  from  Equations  5.12  and 

5.13,  and  calculated  from  Equation  5.30.  For  the  case 

of  the  wing,  the  necessary  parameters  are  given  by 

SS=3  g's  (5.55) 

<*x  = .  1 0  7  6  (5.56) 

px-.2467  (5.57) 

The  results  for  Model  3  are  shown  with  those  of  the 
previous  models  in  Figures  5.8  and  5.9.  The  results  for  the 
different  models  in  both  the  stress  and  range  spaces  are 
much  closer  together  than  the  comparisons  for  the  fuselage 
were.  In  fact,  the  results  of  Figure  5.9  indicate  that  the 
cookie-cutter  approximation  for  wing-loading  is  reasonably 

valid,  especially  for  the  higher  values  of  the  failure 

probability.  There  are  at  least  three  reasons  for  this. 
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First,  the  ratio  of  SK  to  SS  is  aneh  smaller  for  the 


i  U 


i 


wing-loading  case.  The  contiasons  probability  models  rise 
faster  with  increasing  stress,  and  the  differences  between 
the  cookie-cntter  damage  distribution  and  the  continous  ones 
are  decreased.  Second,  the  experimental  results  for  the 
fuselage  showed  the  largest  variation  of  any  of  the  aircraft 
components.  This  makes  the  fuselage  a  worst  case  as  far  as 
reliability  goes.  Finally,  the  aircraft  response  (load 
factor)  is  a  very  steep  function  of  range,  resulting  in  a 
s t e p- f unc t i on  appearance.  All  three  distributions  are  in 
near  agreement  on  SK,  the  sure-kill  range  (point  where 
failure  probability  equals  or  exceeds  .98),  while  the  dif¬ 
ference  in  RS  is  only  about  100  meters. 

Vertical  Stabilizer  Vulnerability.  The  analysis  of  the 
vertical  tail  assembly  can  be  done  in  the  same  way  the  wing¬ 
loading  analysis  was  done.  The  stress  variable  a  is  now  the 
combined  bending  moment  of  the  vertical  tail  and  fuselage. 
The  overpressure  variable  P  is  sampled  as  before,  the  sample 
values  are  propagated  through  the  response  equations,  and 
the  distribution  of  s  is  determined  at  range  points  r. 
Equation  S.19  is  then  integrated  numerically  to  determine 
the  failure  probability  as  a  function  of  range  from  the  1  KT 
standard.  The  stress  model  is  described  as  follows.  For 
every  value  of  P  drawn  at  range  r: 

(1)  Proceed  with  the  response  function  for  wing- 


*.V_V 
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loading  through  step  (5)  to  find  u 


the  component  of  the 


resultant  wind  vector  that  acts  on  the  vertical  tail 


assembly. 

(2)  Determine  the  combined  bending  aoaent  due  to  this 
wind  acting  on  the  vertical  tail  and  fuselage  (a  value  of 


nevton-aeters,  and  SK  specification  as  3.7E6  newton-aeters 


(again,  after  Method  1  of  Reference  [54]).  For  the  cookie- 
cutter  model  of  the  strength  distribution,  one  assuaes  the 
strength  to  be  given  by 

S-b#k  (5.61) 
where  b>k  is  the  sure-kill  value  just  aentioned.  This  being 
the  paraaeter  of  a  Dirac  delta  function,  the  equation  in  the 
range  space  for  the  failure  probability  is  given  by 

pf"1’Fa(bsk)  (5.62) 
where  again  the  distribution  of  the  stress  is  found  using 
nonparaaetr ic  aethods.  The  results  are  displayed  in  Figures 
5.10  and  5.11. 

2«  The  strength  distribution  under  the 
assumptions  of  Model  2  is  a  lognorasl  with  parameters  given 
by  Bridgman's  prescription  [45]  in  Equations  5.7  and  5.8. 
The  SS  and  SK  specifications  just  aentioned  provide  the  S  Q2 
and  S  ob  percentiles  that  go  into  the  calculation. 


The  failure  probability  must  be  found  by  numerically 
integrating  Equation  5.19  at  each  range  point  r.  The 
results  are  displayed  along  vith  those  of  Model  1  in  Figures 
5.10  and  5.11. 

Model  £.  Finally,  under  the  assumptions  of  Model 
3,  the  strength  function  is  given  by  Equation  5.61,  but  the 
failures  are  distributed  as  determined  by  Chenoweth  [20]. 

The  strength  distribution  is  asymptotic  Veibull  vith  para¬ 
meters  given  by  Equations  5.17  and  5.18.  As  usual,  is 

taken  as  1.5  times  the  SS  specification.  From  the  static 
test  results,  ax  and  bx  are  found  to  be: 

ax«1.606  (5.63) 

bx=».143 1  (5.64) 

and  the  SS  specification  is  taken  as  1.763E6  nevton-meters. 

The  results  of  the  analysis  are  displayed  in  Figures 
5.10  and  Figure  5.11  and  are  compared  to  the  other  models. 
Although  there  appears  to  be  some  differences  in  the 
distributions  when  viewed  in  stress  intensity  space,  these 
differences  are  nil  when  viewed  in  range  space.  In 
particular.  Model  3  and  the  cookie-cutter  approximation  are 
indistinguishable.  This  is  a  reasonable  result  since  the 
ratio  of  SK  to  SS  is  low,  and  the  experimental  data  yields  a 
failure  distribution  that  rises  steeply.  From  Figure  5.11, 
RE  is  at  about  330  meters,  while  RS  is  at  about  520  meters. 


The  survivability  of  several  aircraft  components  in  a 
nnclear  blast  environment  has  been  calculated.  Stresa  dis¬ 
tribution  construction  begins  with  stress  functions.  The 
nuclear  weapon  effect  of  intereat  here,  the  overpressure 
variable,  is  a  random  variable  that  enters  into  these  stress 
functions.  The  distribution  of  these  functions,  found  by 
nonpar ame tr ic  methods,  defines  the  stress  distributions  for 
the  components. 

The  strength  distributions  are  surveyed.  Two  common 
approaches,  cookie-cutter  modeling  and  lognormal  modeling, 
are  compared  to  distributions  found  by  analyzing  twenty 
years  of  aircraft  static  test  data. 

The  resulting  failure  probability  with  range  from  the 
weapon  is  calculated  and  compared.  In  general,  one  cannot 
say  beforehand  whether  or  not  a  parametric  model  is  superior 
to  a  cookie-cutter  model.  For  most  cases,  the  results  did 
not  vary  much.  For  the  case  of  the  fuselage,  however, 
neither  the  cookie-cutter  nor  the  lognormal  was  as  conserva¬ 
tive  as  the  actual  test  data  suggested. 

In  the  next  Chapter,  the  survivability  to  nuclear 
thermal  effects  will  be  considered,  and  the  difficulties 
posed  by  the  lack  of  failure  data  will  be  noted. 


VI.  Aircraft  Vulnerability  In  Unclear  Thermal  Environment* 

Overview 

In  this  chapter,  a  s tr e s s- s tr ength  interference  theory 
calculation  of  the  probability  of  failure  of  skin  panels  in 
nuclear  thermal  environments  is  accomplished.  A  search  of 
the  database  for  nuclear  effects  on  aircraft  yields  no 
direct  information  on  either  stress  or  strength  distribu¬ 
tions  for  the  failure  modes  of  interest.  The  temperature  of 
a  thin  skin  panel  in  a  nuclear  thermal  environment  is  a 
fundamental  stress.  Although  this  is  a  deterministic 
response,  statistical  variation  is  introduced  by  treating 
the  radiated  power  from  the  weapon  as  a  random  variable,  and 
then  finding  the  distribution  of  the  functions  of  this 
random  variable.  Two  failure  modes  are  considered.  The 
first  is  based  on  sure-safe  and  sure-kill  specifications  for 
aircraft,  and  is  analogous  to  the  work  done  on  blast  vulner¬ 
ability.  For  this  mode,  two  strength  distributions  are 
examined — a  cookie-cutter  and  a  lognormal.  An  alternate 

failure  mode  is  investigated  that  considers  skin  panel 
yielding  as  a  function  of  both  gust  and  thermal  loading. 
The  calculation  illustrates  the  use  of  the  nonp a r a m e t r i c 
interference  theory  technique  in  treating  problems  where 


multiple  nuclear  effects  contribute  to  a  failure  mode. 


A  Data  Search  Por  Thermal  Effects  on  Aircraft 

As  in  the  case  of  blast  effects  on  aircraft,  no  direct 
failure  data  for  aircraft  components  were  found.  Items  of 
interest  in  the  nuclear  tests  examined  include  thermal  fin¬ 
ance  and  temperature  time  histories.  Peak  temperatures  for 
the  time  histories  of  Operation  &EDVIN6  are  well  below  melt 
point.  A  typical  result  from  Reference  [70]  is  shown  in 
Figure  6.1.  Lacking  actual  service  histories  for  the  fail¬ 
ure  mode  of  interest,  the  only  alternative  is  engineering 
modeling.  An  important  variable,  as  displayed  in  Figure 
6.1,  is  the  temperature  of  a  thin  panel  as  a  function  of 
time.  This  is  a  fundamental  thermal  stress  on  the  system, 
and  is  discussed  in  detail  below. 

A  Deterministic  Thermal  Stress  Model 

A  simple  model  of  the  temperature  rise  in  a  thin-skin 
assembly  can  be  constructed  by  requiring  an  energy  balance 
condition  [43].  That  is,  the  temperature  in  the  thin  skin 
obeys  the  differential  equation: 

pCvAxdT/dt=Yth(t)/ (4nr2)-h(v) (T-Tq)  (6.1) 

where  the  term  on  the  left  of  the  equals  sign  is  the  time 
rate  of  change  of  the  energy  deposition  in  the  material. 
The  variables  in  this  term  include  p,  which  is  the  material 
density  in  kg/m^;  Cy  is  the  specific  heat  in  Joul e s / (kg-°- 
Kelvin);  Ax  is  the  skin  thickness  in  meters)  T  is  the  skin 


temperature  in  degrees  Kelvin;  t  is  the  time  in  seconds. 


The  first  tern  to  the  right  of  the  equals  sign  is  the 
source  term,  corrected  for  spherical  divergence.  The  vari¬ 
ables  there  are  given  by: 

Tth(t)=Pm(W)f <t)Trabco**(t)  (6.2) 
where  PB(W)  is  the  maximum  radiated  power  (in  watts)  from 

the  weapon  and  is  a  function  of  W,  the  weapon  yield  in 
kilotons  (KT);  f(t)  is  the  fractional  radiated  power  at  time 
t;  T£  is  a  transmittance  factor  to  account  for  attenuation 
through  the  atmosphere;  is  the  absorptivity  of  the  skin 
surface;  9  is  the  look  angle  between  the  weapon  and  the 
aircraft  skin  and  is  a  function  of  time  due  to  the  motion  of 
the  aircraft.  The  variable  r  in  Equation  6.1  is  the  slant 
range  (in  meters)  from  the  weapon  at  time  t. 

The  last  term  in  Equation  6.1  is  the  sink  term,  the 
rate  of  energy  dissipated  to  the  environment  by  convective 
cooling.  The  variables  here  include  h(v),  which  is  the 
(velocity  dependent)  heat  transfer  coefficient  in  Joules/(°- 
Kelvin-meter  -sec)  and  T q,  the  temperature  of  the  air  adja¬ 
cent  to  the  skin  in  degrees  Kelvin. 

The  equation  can  be  made  dimensionless  by  introducing 
some  natural  variables.  Since  the  differential  equation  is 
one  that  couples  temperature  and  time,  the  natural  variables 
involve  both  times  and  temperatures.  The  natural  time 
variables  include  tm,  the  time  (in  sec)  to  peak  thermal 
radiated  power.  The  dimensionless  time  t  can  then  be 


defined  as 


Time  ,  seconds 


Tims  v  seconds 

Exposed  facing  tempe.  lure-time  histories  of  NAA  specimen  with  y,-inch 
cell  siae  and  0.016-inch  white-painted  skin.  (Computed  peak  temperature:  243*  F.) 

Figure  6.1.  Temperature  Time  History  In  A  Thin 
Plate  (From  Reference  [70]} 


T  =  t/t 


(6.3) 
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Another  important  time  variable  is  tc,  the  cooling  time, 
given  by 

t0-pCy&*/h.  (6.4) 

where  the  right  hand  side  variables  have  already  been 
described.  One  can  show  that  in  the  absence  of  any  source 
term,  the  cooling  time  is  the  time  for  a  preheated  structure 
to  cool  to  1/e  of  its  initial  temperature,  where  e  is  the 
natural  logarithmic  base.  Finally,  the  ratio  of  time 
constants  will  be  denoted  by 

>  =  tm/tc  (6.5) 

The  natural  temperature  variables  include  the  ambient 
temperature,  Tq ,  in  degrees  Kelvin.  The  dimensionless  time 
dependent  temperature  T  can  be  defined  by 

T(t)=T(t)/T0  (6.6) 

A  time  dependent  characteristic  temperature  can  also  be 
defined.  This  temperature  depends  on  geometrical  and 
physical  characteristics  of  the  skin  material,  and  on 
properties  of  the  thermal  pulse.  The  dimensionless 
characteristic  temperature  can  be  written  as 

Tc(T)=Tk(T) f (t)  (6.7) 

A 

where  Tk(r)  is  given  by 

Tk(T)=Pm(w)Trabcos(9(T))tm/(4',r2pCvAlT0)  (6.8) 

The  variables  on  the  right  hand  side  of  Equation  6.8  have 
been  described  previously.  The  term  f(T)  in  Equation  6.7 
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is  the  relative  radiated  power  at  time  T.  That  is 
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f(r>«P(T>/P_(l>  (6.9) 

m 

where  P  ( r)  is  the  radiated  power  in  watts  at  time  T.  For 

A 

stationary  targets,  T  ^  reduces  to  a  time-independent 
constant.  Dividing  all  times  in  Equation  6.1  by  tm,  the 
time  to  thermal  maximum,  and  all  temperatures  by  Tq,  the 
ambient  temperature,  leads  to  a  dimensionless  differential 
equation.  Solving  that  equation  by  using  an  integrating 
factor  and  the  natural  variables  just  described  leads  to  the 
solution 

T(T)-1  +  [T(0)-1+J^  Tc(r»)e‘yT’dT']exp(-7T)  (6.10) 

Equation  6.10  is  the  deterministic  equation  for  the 
stress  function.  It  represents  the  anticipated  response  of 
the  skin  to  a  nuclear  thermal  pulse.  The  statistical  distri¬ 
bution  of  this  response  can  be  determined  by  propagating 
statistical  information  about  the  input  variables  through 
the  function.  Of  the  many  input  variables  one  might  con¬ 
sider,  the  radiated  power  from  the  weapon  will  be  examined 
statistically  and  the  others  left  as  deterministic.  Two 
reasons  dictate  this  choice.  First,  it  will  provide  a 
parallel  development  to  the  work  of  the  previous  chapter. 
In  particular,  a  thermal  analog  to  the  work  of  Carpenter  and 
Kuhl  [47]  will  be  developed.  Second,  it  is  difficult  to 
find  statistical  data  for  the  other  variables  of  interest. 
In  fact,  the  statistical  information  for  the  radiated  power 
is  not  readily  available,  and  one  must  analyze  the  data  that 


one  can  find.  In  the  next  section,  the  deterministic  stress 


function  is  expanded  into  a  random  variable  equation. 

The  The rma  1  Stress  Distribution  For  A  i  KT  Sea  Level  Bur s t 

General  Approach.  The  objective  now  is  to  find  the 
statistical  description  of  the  temperature  in  a  thin  skin  at 
any  time  desired.  Equation  6.10  is  the  basic  deterministic 
temperature  response.  The  dependence  on  the  basic  weapon 
effect,  the  radiated  power,  can  be  explicitly  demonstrated 
by  writing  the  characteristic  temperature  in  the  form  shown 
in  Equation  6.7.  The  resulting  integral  equation  is: 

T(T)=1+[T(0)-1+J~J"  Tk(T')f(T')e7r,dT']exp(-7r)  (6.11) 

If  one  now  recognizes  that  the  radiated  power,  f(T),  is 

statistically  distributed,  then  the  resultant  output  vari- 

A 

able,  T(t)  is  also  statistically  distributed,  since  it  is  a 
function  of  a  random  variable.  In  particular,  if  the  dis¬ 
tribution  of  the  radiated  power  has  been  determined,  then 

the  particular  sample  value  f^(T)  results  in  a  corresponding 

A 

sample  value  T ±(t)  where  the  two  valuea  are  functionally 
related  through  Equation  6.11.  Given  the  distribution  of 
the  radiated  power,  the  distribution  of  the  temperature  in 
the  structure  can  be  determined  by  the  methods  already 
discussed.  The  determination  of  the  distribution  of  the 
radiated  power  is  the  next  task  to  be  accomplished. 

A  Siaii.jil.iiAl  Mo 4.2.1  01  Jhe,  RadiAiii  Power.  As  previ¬ 
ously  mentioned,  a  search  of  the  database  for  thermal 
effects  on  aircraft  did  not  yield  any  information  about 


mm 
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statistical  thermal  source  models — that  is,  no  analog  to  the 
work  of  Carpenter  and  Kuhl  [47]  was  found.  In  order  to 
continue  with  an  interference  theory  calculation  one  must 
have  a  source  of  data.  There  are  several  approaches  that 
could  be  taken,  but  all  depend  in  part  on  a  deterministic 
model  of  the  weapon  effect.  Two  such  deterministic  models 
are  discussed  below. 

A  relatively  simple  model  of  the  radiated  power 
produced  by  a  nuclear  weapon  is  described  by  Glasstone  [62]. 
Figure  6.2,  reproduced  from  Reference  [62],  displays  the 
normalized  power  as  a  function  of  normalized  time.  In  this 
model,  the  time  to  (second)  thermal  maximum,  and  the  peak 
radiated  power  at  that  time  are  scaled  according  to  yield. 
Consequently,  the  plot  of  Figure  6.2  has  been  used  for  a 
variety  of  yields  and  burst  altitudes.  Glasstone  quotes  the 
results  as  being  valid  +  25%  for  altitudes  below  about  4-5 
kilometers  (km)  and  +  50%  for  altitudes  above  that.  This 
statement  acknowledges  that  the  value  of  the  radiated  power 
is  statistically  uncertain. 

Glasstone's  pulse  is  synthesized  from  both  experimental 
data  gathered  during  the  days  of  atmospheric  testing  and 
theoretical  modeling.  With  the  signing  of  the  Limited  Test 
Ban  Treaty  in  1963,  theoretical  calculations  of  thermal 
environments  have  largely  replaced  data  collection 
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Figure  6.2.  Glasstone's  Thermal  Pulse  (From  Reference  [62]) 


Sharp  of  the  Air  Force  Weapons  Laboratory  (AFWL) 
published  an  alternate  model  in  1973  [77].  The  model  was 
published  as  a  FORTRAN  subroutine,  and  includes  variation  of 
the  thermal  pulse  shape  with  altitude.  This  feature  was  not 
available  in  the  Glasstone  model.  This  AFWL  thermal  power 
routine,  called  FLUX,  is  still  in  use.  It  is  found,  for 
example,  in  some  computer  codes  that  deal  with  nuclear 
survivability  problems.  These  codes  include  TRAP  [65], 
QUANTA  [38],  and  FLEE. 2  [36]. 

The  problem  of  finding  a  statistical  description  of  the 
radiated  power  from  a  nuclear  weapon  still  remains.  One 
approach  which  might  be  used  is  documented  by  Ostermann  and 
Collins  [75].  These  authors  specify  a  lognormal  distribu¬ 
tion  for  a  nuclear  weapon  environment,  using  the  output  of  a 
nuclear  effects  algorithm  as  the  median  of  the  distribution. 
The  scale  parameter  of  the  lognormal  (which  is  just  the 
standard  deviation  of  the  log  of  the  data)  is  said  to  be 
estimated  by  a  combination  of  judgment  and  data.  This 
technique  is  similar  to  Bridgman's  approach  to  specifying 
the  strength  distribution  [45]. 

As  an  alternate  approach,  one  can  examine  the  data  that 
does  exist,  however  sparse,  and  formulate  the  distribution 
of  the  radiated  power  from  that  alone.  Choosing  this  as  a 
rationale,  the  Glasstone  model  cannot  be  used.  Although  an 
uncertainty  is  stated,  it  is  not  clear  how  to  translate  ±.25% 
into  a  statistical  distribution.  The  AFWL  FLUX  routine 


f 


|  statistical  radiated  power  model  based  on  the  data  alone. 

Before  proceeding  with  this,  FLDX  must  be  examined  in  a 
t  little  more  detail. 

I  FLUX  is  a  numerical  model  that  seeks  to  duplicate  the 

thermal  power  output  of  a  more  complex  code  called  SPUTTER 

[80].  SPUTTER  is  a  program  used  to  dynamically  model 

j  nuclear  fireball  phenomenology.  Among  the  quantities 

available  from  a  SPUTTER  run  are  the  thermal  power  and 

energy  as  a  function  of  wavelength  and  time.  Sharp  used 
I 

1  SPUTTER  to  obtain  t i m e -de p e nd e nt  power  and  energy  in  35 

wavelength  bands  covering  the  atmospheric  transmission 
window.  FLUX  was  then  developed  as  a  curve  fitting  tech- 

I  u 

1  nique  to  match  the  unattenuated  source  power  and  energy 

values  as  functions  of  time  summed  over  the  35  wavelength 
bands . 

^  FLUX  therefore  derives  its  sole  credibility  from  the 

SPUTTER  runs  which  it  seeks  to  match.  The  statistical  model 
to  be  developed  is  thus  again  conditional.  SPUTTER  is 
presumed  to  be  true. 

If  the  SPUTTER  runs  are  taken  to  represent  nature,  then 
one  needs  to  assess  statistically  the  ability  of  FLUX  to 
model  SPUTTER  output.  The  answer  depends  on  the  height  of 
burst  and  yield  regime.  Figure  6.3  taken  from  Reference 
[77]  illustrates  a  rather  good  result,  while  Figure  6.4 
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Figure  6.3.  Thermal  Power  Versus  Time  —  198  Kilotons 
At  Sea  Level  (From  Reference  [77]) 
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Figure  6.4.  Thermal  Power  Versus  Time  -  3.8  Megatons  At 
Sea  Level  (From  Reference  [77]) 
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illustrates  a  soaevhat  poorer  one 


Of  special  interest  is 


Figure  6.5,  which  is  the  fit  for  a  1  KT  burst  in  a  sea  level 
atmosphere.  In  these  figures  the  solid  line  is  the  FLUX 
prediction  line,  while  the  circles  are  the  SPUTTER  output. 
The  objective  is  to  find  a  statistical  measure  of  the  qual¬ 
ity  of  the  FLUX  fit. 

One  could  make  the  lognormal  assumption  as  done  by 
Ostermann  and  Collins  [75],  and  use  FLUX  output  as  the 
median  radiated  power  value.  The  scale  parameter  for  the 
lognormal  could  be  estimated  by  the  average  root  mean  square 
error.  That  is,  one  might  take  0p  as: 

Pp“Si“l(lnIf-lnyi)2/n}1/2  (6.12) 

where  y^  is  a  particular  SPUTTER  value  and  Yf  is  the  FLUX 
prediction  at  the  same  time.  The  only  problem  with  this  is 
that  there  is  no  guarantee  that  the  FLUX  line  is  the  median 
line.  This  is  borne  out  by  a  second  examination  of  Figure 
6.5.  One  sees  that  FLUX  und  e  r  pr  e  d  i  c  t  s  SPUTTER  in  the  late 
time  regime  in  particular.  A  calculation  of  the  mean 
residual,  1,  defined  by: 

S2i  =  l(lnYf“layi)/n  (6.13) 
for  the  FLUX  fit  to  SPUTTER  run  FB-21  (Figure  6.5)  yields  a 
value  of  -.18,  showing  the  overall  und e r p r e d i c t i on.  The 
scale  parameter  as  estimated  in  Equation  6.12  is  about  .50. 
This  fit  can  be  improved  in  the  following  way. 


J. 
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Since  the  objective  of  FLUX  is  to  sitch  SPUTTER  point 
for  point,  FLUX  may  be  thought  of  as  a  predictor  variable 
for  the  SPUTTER  output.  Hence,  if  FLUX  exactly  matches 
SPUTTER,  then  a  plot  of  SPUTTER  versus  FLUX  for  fixed  yield 
and  height  of  burst  will  yield  a  straight  line  vith  unit 
slope  and  zero  intercept.  The  SPUTTER  data  for  the  1  XT  sea 
level  case,  vith  the  corresponding  FLUX  prediction  values  is 
shown  in  Table  VI. 1.  In  this  table,  time  enters  in  as  an 
implicit  variable.  The  regression  of  SPUTTER  on  FLUX  is 
shown  in  Figure  6.6.  As  already  noted,  at  late  times  FLUX 
underestimates  the  SPUTTER  values.  This  is  seen  in  Figure 
6.6  in  the  low  power  regime.  The  one  sigma  boundaries  are 
also  shown  on  the  plot.  Linearity  is  pronounced,  as 
expected,  with  a  correlation  coefficient  of  .96.  The 
regression  equations  in  the  log  space  are  given  by: 

<lnYs>=mlnYf+b  (6.14) 

ffS|F  =  *3258  (6.15) 

where  Ys  is  the  median  SPUTTER  output  in  watts,  Y£  is  the 
FLUX  prediction  in  watts,  m  is  .7626,  b  is  6.847,  and  ^|p 
is  the  conditional  variance  of  the  SPUTTER  output  given  FLUX 
as  a  predictor  variable.  As  a  comparison,  the  mean  residual 
for  this  fit  has  been  reduced  to  about  8E-6,  while  the  root 
mean  square  error  has  been  reduced  to  about  .30.  The  least 
squares  fit  incorporates  the  standard  assumption  that  the 
distribution  of  the  residuals  is  normal. 
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TABLE  VI. 1 


DATA  FOE  1  KILOTON  AT  SEA  LEVEL 


Power 

(Watts) 


FLUX  Prediction 
(Watts) 


1.0 

2 . 444E11 

8 .668E10 

0.9 

2 . 977E11 

1 .043E11 

0 . 8 

3 . 660E11 

1 .284E11 

0.7 

4.539E11 

1 .624E11 

0.6 

5 . 5  8 1 El 1 

2 . 130E11 

0.5 

6 .709E11 

2.93  5E1 1 

0.4 

8.049E11 

4 . 345E11 

0.3 

1 .000E12 

7 . 201E1 1 

0.2 

1 .383E12 

1 .462E12 

0.15 

1 .792E12 

2.3  96E1 2 

0.1 

2 . 744E12 

4 . 604E1 2 

0.09 

3 . 1 6  6  El  2 

5 . 343E12 

0.08 

3 . 864E12 

6 .193E12 

0.07 

5 .056E12 

7 . 13  4E12 

0.06 

8 . 269E12 

8.594E12 

0.05 

1 .223E13 

1 .200E13 

0.04 

1.402E13 

1 .466E13 

0.03 

1.348E13 

1.374E13 

0.02 

1 .019E13 

8 . 504E1 2 

0.015 

6 . 941E12 

5  . 10  5E12 

0.01 

3  . 392E12 

2 . 3  87E12 

0.009 

2 . 804E12 

1 .992E12 

0.008 

2 . 342E12 

1 .651E12 

0.007 

1.909E12 

1 .364E12 

0.006 

1.525E12 

1 .132E12 

0.005 

1 .185E12 

9 .539E11 

0.004 

8.886E11 

8 .347E11 

0.003 

7 . 248E11 

7 . 960E11 

0.0026 

7 . 07  3E11 

8 . 188E11 

0.002 

8.242E11 

9 . 129E11 

0.0015 

9 . 794E11 

1 .092E12 

0.001 

9.851E11 

1 .602E12 
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SPUTTER  on  Ln  FLUX 


Since  the  plot  is  in  the  log  space,  the  radiated  power 
Ig  is  a  lognormal  random  variable  with  location  parameter  ap 
and  scale  parameter  0p  given  by: 

ap=<lnYs>  (6.16) 

Pp=<TS|F  (6*17) 

The  statistical  model  of  the  radiated  power  is  thus 

determined.  The  explicit  expression  for  a  sample  value  of 
the  relative  radiated  power  can  be  written  in  terms  of  the 
above  parameters.  Given  z^,  a  value  of  the  standard  normal 
variate,  fi(T)  is  given  by: 

fi(T)=eip{<»p(r)+PpZi)  (6.18) 
The  i  subscript  denotes  the  statistical  meaning  of  the 
particular  value  of  f^  in  terms  of  the  standard  normal 
variate  z^.  In  particular, 

Pr  (Fif  ■  }=Pr(ZiziJ  (6.19) 
Equation  6.18  may  be  directly  inserted  into  Equation  6.11  to 
complete  the  expression  for  the  distribution  of  the 
temperature  in  a  thin  skin  assembly. 

With  a  statistical  model  for  the  temperature  in  a  thin 
skin  completed,  the  stress  distribution  is  known.  The  other 
half  of  the  problem  is  the  determination  of  the  strength 
distribution.  The  choice  of  the  strength  variable  depends 
on  the  failure  mode  of  interest.  Melt  mode  failures  will  be 


examined  in  the  next  section 


The rma 1  Vulne  r ab i 1 i t v  Model in 


The  statistical  description  of  the  temperature  in  a 
thin  shin  has  been  completed.  This  may  be  the  stress  vari¬ 
able  of  interest,  bnt  it  may  also  appear  in  the  strength 
function,  depending  on  the  failure  model  chosen.  The  choice 
of  the  strength  variable  also  depends  on  the  failure  mode  of 
interest.  Two  failure  models  are  considered  below.  In  one, 
melt  mode  failures  are  considered.  In  the  other,  the  tem¬ 
perature  of  the  shin  enters  into  both  the  stress  and 
strength  distributions. 

Me.iJt-Mode.  Vulnerability.  In  the  previous  chapter, 
methods  from  the  nuclear  survivability  literature  [54,45] 
were  compared  with  data  that  allowed  one  to  infer  actual 
strength  distributions  for  several  aircraft  components. 
Aircraft  must  be  designed  for  a  certain  amount  of  gust 
hardness  since  atmospheric  turbulence  is  a  natural  aircraft 
environment.  It  is  not  surprising,  then,  that  enough  data 
exists  in  regards  to  mechanical  loading  to  formulate 
strength  distributions  directly  from  test  data  [22].  Some 
have  even  suggested  that  because  of  this,  modern  day  air¬ 
craft  have  some  inherent  hardness  to  nuclear  gusts  [72]. 

In  contrast,  aircraft  are  not  normally  designed  for 
thermal  environments  lihe  that  from  a  nuclear  weapon.  Con¬ 
sequently,  no  direct  data  exists  from  which  to  infer  a 
strength  distribution.  There  is  then  no  choice  but  to  model 
the  strength  function  by  choosing  values  of  skin  tempera- 
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tares  that  are  considered  critical.  In  Reference  [54]  two 
such  values  are  chosen — a  sure-safe  (SS)  temperature,  Tgg, 
and  a  sure-kill  (SK)  temperature,  Tgg*  The  sure-safe 
temperature  is  that  temperature  which  results  in  a  20% 
reduction  in  the  modulus  of  elasticity.  The  sure-kill  tem¬ 
perature  is  the  melt-point  of  the  material.  Two  strength 
distributions  will  be  considered.  As  in  the  case  for  blast, 
a  cookie-cutter  strength  distribution  will  be  called  Model  1 
while  the  a  priori  lognormal  strength  distribution  will  be 
referred  to  as  Model  2. 

Cookie-Cutter  Failure  Distribution.  The  structure 
fails  upon  encountering  the  SK  temperature,  and  does  not 
fail  otherwise — i.e,  the  failure  probability  density  func¬ 
tion  (PDF)  is  a  Dirac  delta  function  centered  at  Tgg. 

A  Priori  Log.no  rma.i  Failure  Digirib  uiio  n .  The 
probability  of  failure  for  the  skin  structure  is  .98  if  the 
SK  temperature  is  encountered;  it  is  .02  if  the  SS  stress 
value  is  encountered;  the  distribution  of  failures  is  log¬ 
normal  [45]. 

These  two  distributions  are  illustrated  in  Figure  6.7. 
The  stress  space  is  now  temperature.  The  circles  represent 
Model  2,  while  the  solid  line  represents  Model  1.  Model  1 
yields  the  familiar  cookie-cutter  plot  in  the  absence  of  any 
statistical  model  of  the  stress.  However,  since  a  statisti¬ 
cal  model  of  the  stress  does  exist,  the  failure  probability 
will  be  a  continuous  function  of  range.  For  Model  1,  it  is 


given  by: 

Pj(r)ml-Fj(Tg£)  (6.20) 

where  Fg(T)  is  the  cumulative  distribution  function  (CDF) 
of  the  temperature  (the  stress  variable)  in  the  thin  skin  at 
the  time  of  interest. 

The  resultant  vulnerability  of  the  skin  is  shown  in 
Figure  6.8  for  both  Model  1  and  Model  2  strength  distribu¬ 
tions.  The  c o o k i e - c u t t e r  damage  distribution  is  not  very 
different  from  the  lognormal  one.  The  failure  probability 
is  a  steep  function  of  range,  dropping  from  .98  to  .02  in 
the  span  of  100  meters  even  for  the  lognormal  assumption. 
The  reason  for  this  is  the  rapid  decrease  in  peak  tempera¬ 
ture  with  range.  This  is  true  even  if  the  structure  is  not 
cooled.  If  it  is,  the  decrease  is  even  faster.  Thus,  one 
sees  that  cookie-cut  ter  techaiqnes  can  work  well  depending 
on  the  rate  of  change  of  the  critical  response  variable  with 
other  variables  in  the  problem. 

Before  leaving  the  thermal  vulnerability  problem,  one 
can  consider  an  alternate  failure  mode.  In  particular,  one 
can  examine  the  effects  of  combined  gust/thermal  loading  on 
the  mechanical  integrity  of  a  structure.  This  is  the  last 
topic,  and  is  considered  in  the  next  section. 

Comb ine  d  Blast/ The  r ma 1  Vulnerab il itv.  The  fundamental 
idea  behind  most  thermal  vulnerability  calculations  is  to 
find  the  amount  of  heat  required  for  some  specific  effect. 
The  topic  of  yielding  of  skin  panels  due  to  thermal  loads 
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combined  with  internal  pressurization  loads  is  treated  using 
Method  2  of  Reference  [54].  This  problem  demonstrates  that 
a  combination  of  conditions  can  precipitate  failure  rather 
than  the  primary  weapon  effects  acting  alone.  On  the  other 
hand,  a  given  failure  mode  might  depend  on  more  than  one 
weapon  effect.  In  the  alternate  failure  model  described 
below,  the  combined  effects  of  blast  and  thermal  in  causing 
yielding  of  wing  skin  panels  is  investigated.  This  will 
illustrate  the  utility  of  the  stress-strength  interference 
theory  technique  in  treating  such  problems. 

S_lX.es.  A.  On  A  Skin  Pa.ne.1..  A  stress  analyst  would 
attempt  to  calculate  the  stresses  in  the  skin  rather  than 
use  the  SS  and  SK  specifications  of  the  previous  section 
[43].  A  representative  wing  section  is  shown  in  Figure  6.9. 
Treated  as  a  simple  box  beam  with  constant  bending  stress 
[76],  and  using  a  thermal  model  from  a  classical  text  like 
Gatewood  [61],  one  can  show  that  the  stress  in  the  lower 
section  at  position  z  and  time  t  is  given  by: 

Ox(x,t)-flrg(z,t)+ath(z,t)  (6.21) 
where  <yg(z,t)  is  the  stress  component  due  to  the  gust  load¬ 
ing  and  ff^Cz.t)  is  the  stress  component  due  to  the  thermal 
loading  from  the  nuclear  detonation.  These  two  terms  are 
given  by: 

<jrg(z,t)  —  <*lgNLz/c  (6.22) 

^th^,t)  =  <^thl(z,t)*^th2^z’^^  ^t  h3^z’^  (6.23) 
The  thermal  stress  terms  [61]  are  given  by: 
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(6.24) 


ffthl(z*  t)“_alE(T>AT(z't) 

*th2(z' t)=(albwE(T)/Ax)/-c-e4T(z*t)dz  (6.2  5) 

tfth3(z* t)  =  (olbwzE(T)/Ix)/-c-tzAT(z' t)dz  (6.2  6) 

The  variable*  in  the  above  equations  are  described  as 
follows.  g  is  the  allowable  design  stress,  taken  as 

6.891E7  pascals  (10*  psi).  is  the  load  factor  on  the 

wings  (dimensionless).  The  coordinate  z  is  the  location 
where  the  stress  is  desired  (meters)  measured  relative  to 
the  center  of  the  assembly.  The  variable  c  is  the  midpoint 
of  the  skin  as  measured  from  the  center  and  is  taken  as 
.1524  meters  (6  inches).  6  is  the  h a  1 f- 1 h i c kne s s  of  the 
skin,  taken  as  1.588E-3  meters  (.0625  inch).  a j  is  the 

coefficient  of  linear  expansion  (  2.286E-5  ■/■-“-Kelvin). 
E(T)  is  the  modulus  of  elasticity  (in  pascals)  at 
temperature  T  (in  0  -Kelvin).  The  room  temperature  modulus 

7 

of  aluminum  is  taken  as  6.891E10  pascals  (10  psi).  As 
previously  discussed,  T  is  the  temperature  of  the  skin  in 
degrees  Kelvin.  The  term  AT  is  given  by: 

AT(z,t)*T(z,t)-T0(z.t)  (6.27) 

where  T(z,t)  is  the  temperature  at  location  z  and  time  t  and 
Tq  is  the  ambient  temperature  in  degrees  Kelvin.  The  dimen¬ 
sion  bw,  as  shown  in  the  figure,  is  the  chordwise  length, 

taken  as  .9144  meter  (36  inches).  Finally,  Ax  is  the  cross- 

2 

sectional  area  (15  square  inches  or  9.677E-3  meter  ),  while 
Iz  is  the  area  moment  of  inertia  (474  in*  or  1.973E-4 
meter*)  . 
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Nov,  using  the  thin  skin  approximation  as  before,  the 
temperature  response  of  the  skin  is  as  shovn  in  the  previous 
section.  The  bottom  skin  panel  is  heated  uniformly  to  a 
peak  temperature  Tp  so  that  the  temperature  profile  with  z 
is  given  by: 

T(z)=Tp  v  -c-Mzl-c+t  (6.28) 

T(z)=Tq  elsewhere  (6.29) 

where  t  is  the  skin  hal f- thi ckne s s ,  as  previously  defined. 

With  the  above  temperature  profile,  the  integrals  in 
Equations  6.25  and  6.26  are  calculable  in  closed  form.  The 
results  for  the  lower  wing  station  (z--c)  for  the  geometry 
of  Figure  6.9  lead  to  the  equation: 

ffx=<rlgNL-.3582a1E(Tp)ATp  (6.30) 

Equation  6.30  is  a  function  of  two  random  variables-- 
N^,  the  load  factor  on  the  wings,  and  Tp,  the  peak  tempera¬ 
ture  in  the  thin  skin.  This  equation  is  the  stress  response 
of  the  skin.  The  response  is  dependent  on  both  the  gust 
loading  and  the  thermal  loading  caused  by  the  nuclear  wea¬ 
pon.  Even  though  this  is  a  function  of  two  random  variables, 
it  is  not  clear  that  it  needs  to  be  treated  that  way.  The 
gust  loading  and  peak  thermal  loading  can  appear  at  very 
different  times,  and  this  would  effectively  decouple  the  two 
effects . 

The  relative  vulnerabilities  of  the  four  parts  of  the 
box  beam  are  discussed  in  more  detail  in  Appendix  E.  In 
fact,  the  analysis  there  shows  that  the  lower  skin  is  not 
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the  most  vulnerable  assembly.  However,  that  part  of  the 


beam  is  interesting  in  that  the  gust  dominates  the  stress 
distribution,  while  the  thermal  dominates  the  strength 
distribution.  Since  the  thermal  contribution  to  the  load 
stress  is  slight,  one  can  first  consider  the  case  of  gust 
loading  only. 

The.  Case  of  Gust  Loading  On l.y .  As  a  first 
approximation.  Equation  6.30  is  evaluated  when  the  thermal 
environment  is  absent.  This  reduces  the  stress  function  to: 


Wl  <6-31> 

The  strength  function  is  taken  as  the  yield  point  of 


6061-T6  aluminum.  That  is,  the  critical,  or  sure-kill 
stress  value  is  given  by: 

<Tsk=2.756E8  pascals  (40,000  psi)  (6.32) 

Since  this  is  a  c o o k i e - c u t t e r  strength  distribution, 
the  reliability  interference  integral  reduce^  to  the  result: 


Pf  =  l-Fm(<Tsk)  (6.33) 
where  a  is  the  stress  variable,  and  is  just  <rx  as  given  in 
Equation  6.31.  The  overlay  of  the  wing  stress  model  with 
the  previously  considered  cookie-cutter  and  lognormal  gust 
vulnerability  models  are  shown  in  Figure  6.10.  Model  1  is 
the  cookie-cutter  model  for  gust  loads  discussed  in  Chapter 
V,  while  Model  2  is  the  a  priori  lognormal  model,  and  Model 
3  is  the  results  based  on  Chenoweth's  data.  The  skin  stress 
model  just  developed  (represented  by  the  Q's  in  the  figure) 
is  a  bit  more  conservative  than  the  other  models  considered. 
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This  is  in  part  because  it  is  j  at  t  another  c o ok i e - c a 1 1 e r 


damage  model,  bat  with  a  lower  sure-kill  specification. 
That  is,  the  snre-hill  load  factor  is  taken  as  4  instead  of 
6.75.  This  value  of  4  is  closer  to  what  the  authors  of 
STEAT- SURVIVOR  [42]  would  call  the  unbiased  cookie-cutter 
damage  specification.  In  any  case  one  sees  that  a  reason¬ 
able  independently  derived  failure  model  gives  about  the 
same  results  as  the  others  considered. 

The  Theraul  S_t£e.s..s  Contribution.  Equation 
6.30  can  also  be  evaluated  with  the  peak  thermal  environment 
present,  and  the  aircraft  in  a  straight  and  level  flight 
condition.  In  that  case,  and  for  the  geometry  previously 
indicated,  the  response  equation  reduces  to: 

<Tx  =  <rlg-.3582a1E(Tp)ATp  (6.34) 

The  thermal  stress  contribution  adds  a  compressive 
stress  term  to  an  existing  tensile  load.  The  thermal  stress 
contribution  is  quite  small  for  the  lower  skin  assembly. 
This  is  because  not  much  heating  has  taken  place  by  the  time 
the  blast  wave  arrives.  For  example,  the  statistical  dis¬ 
tribution  of  temperature  at  blast  arrival  time  (.24  sec)  is 
shown  in  Figure  6.11  for  a  target  at  200  meters  range.  This 
variable  is  statistically  distributed  at  each  time  step  due 
to  the  statistical  variation  in  the  calculation  of  the 
thermal  power  (Figures  6. 3-6. 3).  The  position  at  200  meters 
range  is  well  inside  the  sure-kill  region  for  gust  acting 
alone.  The  illustration  in  Figure  6.11  shows  a  temperature 
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maximum  in  the  right  hand  tail  of  the  di  stribat  ion  at  less 
than  585  °. 

This  does  not  mean  that  a  gust/thermal  synergism  does 
not  exist.  It  only  means  that  it  does  not  exist  in  the 


stress  function  for  the  lower  skin.  The  behavior  of  the 

strength  function  there  is  examined  next. 

Gas  t  Dependent  Stress-Thermal  Dependent  Strength. 

Even  though  the  thermal  loading  does  not  contribute  much  to 

the  stress  function,  a  combined  effects  problem  can  still 

exist.  The  yield  point  of  the  material  has  been  treated  so 

far  as  a  constant.  Hence,  the  strength  distribution  has 

been  considered  as  a  c o o k i e - c u t t e r  in  the  stress  space. 

However,  the  yield  point  can  vary  with  temperature.  A 

possible  variation  of  yield  strength  with  temperature  is 

illustrated  in  Figure  6.12.  The  authors  of  Reference  [79] 

state  that  the  degradation  of  yield  strength  parallels  that 

o 

of  the  elastic  modulus  out  to  about  480  Kelvin.  At  that 
point,  the  length  of  time  at  a  given  temperature  ("soak 
time")  becomes  an  important  variable  as  well.  Since  480°  is 
a  rarely  observed  temperature  for  the  1  KT  scenario,  the 
strength  function  and  the  elastic  modulus  can  both  be  mod¬ 
eled  by: 

E(T)=E(T0)g(T)  (6.35) 

aSK(T)  =  <TSK(T0>8(T>  (6*36) 

where  g(T)  is  the  approximately  linear  function  shown  by  the 

dashed  line  in  Figure  6.12,  and  is  given  by  Equation 
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When  the  temperature  dependence  for  the  strength  vari¬ 
able  is  taken  into  account.  the  previous  cookie-cutter 
strength  distribution  is  altered.  The  statistical  uncer¬ 
tainty  in  the  temperature  of  the  skin  propagates  into  the 
temperature  dependence  of  the  yield  point  of  the  material.  A 
set  of  calculations  was  performed  to  determine  the  strength 
distributions  from  200  to  500  meters  at  blast  arrival  time. 
The  resulting  continuous  strength  distributions  are  shown  in 
Figures  6.13,  14,  and  15  for  range  positions  of  200,  400, 
and  500  meters  respectively.  The  strength  of  the  part  is 
statistically  distributed  because  the  temperature  of  the 
part  is  (Figure  6.11)  and  the  strength  depends  on  the  tem¬ 
perature  as  shown  in  Equation  6.36.  Since  the  distribution 
of  temperatures  changes  at  each  range  step  the  strength 
distribution  does  also.  For  the  close-in  case  (200  meters) 
the  original  step  function  damage  distribution  of  Equation 
6.32  centered  at  2.756E8  pascals  has  been  transformed  into  a 
continuous  function  with  a  long  left-ward  tail.  The  mid¬ 
range  strength  distribution  (400  meters)  shows  a  cookie- 
cutter  shape  beginning  to  form,  with  the  most  probable 
failure  point  being  the  room  temperature  failure  point  of 
2.756E8  pascals.  The  c o o k  i  e - c u 1 1 e r  limit  has  essentially 
been  reached  in  Figure  6.15  when  the  aircraft  is  500  meters 
from  ground  zero  at  burst  time.  Hence,  even  an  original 
c ook i e - c u t t e r  assumption  (as  in  Equation  6.32)  can  lead  to 
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Figure  6.12.  Degradation  of  Strength  Properties 
With  Instantaneous  Heating  (From  Reference  [79]) 
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continuous  probability  distributions  when  the  statistical 
nature  of  other  variables  is  introduced. 

The  failure  probability  at  each  range  point  depends  on 
more  than  the  strength  distribution  alone.  It  is  the  combi¬ 
nation  of  stress  and  strength  at  every  instant  of  time  that 
governs  the  failure  probability.  The  stress  function  is 
given  by  Equation  6.31  for  the  lover  skin  since  the  gust 
loading  dominates  the  stress,  as  just  discussed.  The  load 
factor,  N^,  is  the  random  variable  that  gives  the  stress 
function  its  random  character.  The  stress  distributions  for 
range  positions  of  200,  400,  and  S00  meters  are  shown 
respectively  in  Figures  6.16  through  6.18. 

With  both  distributions  determined,  one  from  the 
nuclear  blast  effects,  and  the  other  from  the  thermal 
effects,  the  failure  probability  can  be  directly  calculated 
by  the  methods  previously  described  and  illustrated.  The 
failure  probability  as  a  function  of  range  is  displayed  in 
Figure  6.19.  The  results  are  nearly  indistinguishable  from 
the  approximation  that  considered  gust  effects  only  and  the 
room  temperature  value  of  the  yield  stress.  This  is  because 
very  little  heating  has  taken  place  at  gust  arrival  time. 
One  should  not  infer  that  gust/thermal  combinations  are 
never  important,  since  results  may  be  different  for  differ¬ 
ent  parts  of  the  beam  and  for  other  scenarios.  However,  the 
results  do  show  the  utility  of  the  new  nonparame  tr  ic 
approach  to  stress-strength  interference  theory  when  direct 
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Sniaarv 

Aircraft  vulnerability  to  the  nuclear  theraal  pulse 
from  a  low  yield  weapon  has  been  considered.  A  search  of 
the  available  data  sources  for  thermal  effects  on  aircraft 
yielded  no  direct  information  on  either  stress  or  strength 
distributions.  A  deterministic  model  of  the  stress  was 
developed  instead.  The  statistical  variation  in  one  of  the 
input  quantities  to  this  model  was  determined  from  the 
nuclear  effects  literature.  Stress  distributions  were  then 
inferred  by  finding  the  distributions  of  functions  of  this 
random  input  variable.  Melt  mode  failures  were  examined  and 
the  failure  probability  with  range  calculated.  The  results 
show  little  variation  between  a  c o ok i e -cu t t e r  and  a  log¬ 
normal  strength  distribution  model.  A  combined  effects 
problem  was  also  analyzed.  In  this  case  the  blast  environ¬ 
ment  was  found  to  dominate  the  stress  distribution,  while 
the  thermal  effect  determined  the  strength  distribution.  The 
failure  probability  as  a  function  of  range  was  again  calcu¬ 
lated,  with  little  synergistic  effect  noted.  These  problems 
further  illustrate  the  utility  of  the  nonparame tr ic  inter¬ 
ference  theory  technique  in  assessing  nuclear  survivability. 


s 
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VII.  Summary  And  Recommendations 


Summary 

A  new  approach  for  assessing  the  survivability  of 
aircraft  components  in  nnclear  blast  and  thermal 
environments  has  been  presented.  A  nonpar ame tr ic  technique 
for  finding  the  distribution  of  functions  of  random 
variables  has  been  discovered  and  presented.  This  procedure 
allows  one  to  rigorously  determine  the  strength  and  stress 
distributions  for  aircraft  components  exposed  to  nuclear 
blast  and  thermal  environments.  If  direct  service  histories 
are  not  available,  strength  distributions  may  still  be 
inferred  by  considering  the  statistical  variation  in  the 
inputs  to  a  strength  function.  Similarly,  stress 
distributions  may  be  inferred  even  though  direct  stress 
measurements  are  not  possible.  The  reliability  interference 
integral  can  then  be  solved  resulting  in  continuous  failure 
probabilities  as  a  function  of  range  from  a  nuclear  weapon. 
In  the  following  paragraphs,  each  chapter  is  briefly 
reviewed,  and  recommendations  for  future  work  presented. 

Chapter  II  provided  a  review  of  nuclear 
survivabil ity/vulnerabil  ity  methods,  both  deterministic  and 
probabilistic.  Deterministic  methods  do  not  provide  for 
continuous  damage  probabilities  with  range.  Probabilistic 
methods  of  two  broad  types  have  been  used.  Users  of  the 
first  type  model  failure  probabilities  directly  with  range. 
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while  users  of  the  second  model  strength  and  stress 


functions  as  random  variable  processes.  The  latter  is  the 
basic  approach  of  mathematical  reliability  theory. 

Chapter  III  provided  a  brief  review  of  mathematical 
reliability  theory,  including  stress-strength  interference 
theory  and  the  interference  integral.  Engineering  determin¬ 
ism  was  found  to  be  a  special  case  of  this  theory.  That  is. 
c o ok i e - c ut t e r  failure  distributions  result  if  Dirac  delta 
functions  are  used  to  represent  the  probability  density 
functions  for  the  stress  and  strength  random  variables. 
Even  though  these  theoretical  considerations  are  attractive, 
stress-strength  interference  theory  is  difficult  to  apply  to 
large  engineering  systems.  Three  problems  exist — (a)  the 
difficulty  of  developing  a  system  reliability  model  from 
component  models,  (b)  the  analytic  difficulty  of  finding 
distributions  of  functions  of  random  variables,  and  (c)  the 
limited  amount  of  data  available.  Some  of  the  current  meth¬ 
ods  of  approaching  these  difficulties  were  reviewed.  Fault 
tree  analysis,  expectation  analysis,  direct  Monte  Carlo 
simulation,  indirect  Monte  Carlo  simulation,  variable  trans¬ 
formation  techniques,  and  Bayesian  inference  have  all  been 
used  in  attacking  these  problems. 

In  Chapter  IV,  a  new  application  of  nonpa r am e t r i c 
statistics  was  presented  that  allows  one  to  find  the 
distribution  of  functions  of  multiple  random  variables.  The 
technique  can  be  applied  without  using  random  number 
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generators.  Although  nothing  can  be  done  when  little  data 
is  available,  the  nonparaaetric  method  does  provide 
protection  against  drawing  unwarranted  conclusions  from 
sparse  data  bases.  In  addition  to  this,  the  nonp a r a m e t r i c 
approach  eliminates  the  need  for  density  function 
identification,  parameter  estimation,  and  the  taking  of 
partial  derivatives.  An  example  from  the  reactor  safety 
literature  was  presented  to  illustrate  the  method,  and 
provide  a  benchmark  calculation. 

In  Chapter  V,  this  new  theory  was  applied  to  the  prob¬ 
lem  of  aircraft  survivability  in  nuclear  blast  environments. 
In  this  case,  the  strength  distributions  for  several  air¬ 
craft  piece-parts  were  taken  from  the  literature.  A  statis¬ 
tical  model  of  the  overpressure  from  a  nuclear  weapon  was 
also  available  from  the  literature.  This  information 
allowed  the  stress  distribution  to  be  determined  using  the 
new  no np a r a m e t r i c  tool.  The  stress-strength  interference 
integral  was  then  solved,  resulting  in  a  continuous 
probability  of  failure  with  range.  Except  for  the  fuselage, 
cookie-cutter  approximations  appeared  to  be  adequate  for  the 
low-yield  scenario  chosen. 

In  Chapter  VI,  a  more  difficult  problem  was  approached. 
The  analysis  of  the  thermal  vulnerability  of  aircraft  is 
difficult  owing  to  the  lack  of  available  data.  Strength 
distributions  cannot  be  determined  based  on  servioe  histor¬ 
ies.  Stress  distributions  cannot  be  determined  directly 
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either.  However,  stress  and  strength  functions  for  failure 
modes  of  interest  can  still  be  postulated  exactly  as  is  done 
in  the  deterministic  case.  Any  statistical  information  that 
affects  the  stress  and  strength  functions  can  be  incorpo¬ 
rated  using  the  no np a r a m e t r  i  c  technique.  Stress  and 
strength  distributions  are  thus  rigorously  inferred,  based 
on  any  available  data.  The  scatter  in  the  prediction  of  the 
thermal  radiated  power  was  analyzed  in  this  fashion,  and  the 
probability  of  failure  for  shin  panels  examined.  In  this 
case,  cookie-cutter  techniques  were  found  to  be  completely 
adequate  for  the  yield  scenario  chosen. 

Re  commend at  ions  for  Future  fork 

Deterministic  modeling  will  always  be  the  mainstay  of 
nuclear  survivability  assessment.  The  physics  of  nuclear 
weapon  effects  and  the  response  of  weapons  systems  to  those 
effects  will  be  a  topic  of  study  for  years  to  come. 
Probabilistic  modeling  can  and  should  augment  this  work.  If 
data  exists,  however  sparse,  the  n o np a r a m e t r  i  c  tool 
described  in  this  dissertation  can  be  used.  Finding  data 
and  rigorously  processing  it  is  hard  work.  The 
survivability  analyst  must  decide  whether  this  task  is  worth 
it.  Methods  need  to  be  developed  to  assist  in  answering  the 
question,  "Is  this  worth  doing?  Vill  a  'cookie-cutter' 
approximation  be  a  good  one?" 
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One  of  the  more  difficult  problems  in  stress-strength 
interference  theory  is  in  the  assessment  of  very  low  failure 
pr obab i 1 i t e s .  The  nonparame tr ic  tool  applied  in  this 
dissertation  conld  in  principle  be  applied  to  these  problems 
also.  The  assessment  of  a  very  low  failnre  probability 
would  require  a  large  number  of  points  and  substantial 
computer  resources.  Research  in  this  area  should  be 
conducted . 

Vork  should  definitely  continue  in  the  area  of  applied 
n o np a r a m e t r i c  statistics.  This  work  should  involve 
improvements  in  nonp a r a m e t r  i  c  tools  themselves,  and  new 
applications  of  those  tools. 

As  far  as  work  on  the  tools  'themselves,  several 
improvements  need  to  be  made.  Survivability  assessment 
based  on  direct  probability  density  function  (PDF) 
estimation  might  be  useful.  Since  integration  is  well- 
posed,  such  a  tool  would  not  suffer  from  some  of  the  defects 
that  can  arise  when  estimating  a  PDF  by  differentiation  of  a 
cumulative  distribution  function  (CDF).  Endpoint 
extrapolation  is  another  area  that  could  be  improved.  The 
currently  used  center-difference  scheme  for  integrating  the 
PDF  might  be  enhanced  by  using  higher-order  polynomials  or 
spline  methods.  Another  area  of  study  is  stylized  sampling 
from  non-monotonic  functions.  A  preliminary  numerical  inves¬ 
tigation  seemed  to  indicate  that  CDF  estimation  by  stylized 
sampling  worked  well  enough  if  one  was  willing  to  put  up 
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with  jumps  in  the  CDF.  The  size  of  the  jumps  cm  be  made  as 
small  as  desired  by  simply  taking  more  and  more  points. 
However,  the  PDF  estimator  derived  by  differentiating  such  a 
CDF  can  be  a  numerical  problem.  Beyond  this,  a  good 
theoretical  proof  here  would  be  welcome. 

Nonp a r am e t r i c  estimation  techniques  should  also  be 
applied  to  other  problems  in  engineering  and  physics.  One 
simple  yet  potentially  useful  application  might  be  in  Monte 
Carlo  radiation  transport.  Even  if  used  only  as  a  tool  to 
provide  fast  inversion  of  distribution  functions,  computer 
time  might  be  substantially  decreased  in  some  of  these  large 
codes.  Other  applications  will  no  doubt  be  found. 

Finally,  as  applied  to  future  work  in  nuclear  surviva¬ 
bility,  the  most  useful  efforts  would  be  in  determining  the 
statistical  uncertainty  in  nuclear  effects  predictions. 
This  is  a  difficult  task,  involving  a  careful  search  of  all 
known  data;  however,  as  noted  earlier,  if  the  environments 
of  the  nuclear  effects  can  be  statistically  described,  con¬ 
tinuous  failure  probabilities  with  range  will  result,  even 
if  strength  distributions  must  be  taken  as  cookie-cutter. 
Furthermore,  these  environmental  inputs  would  remain  the 
same,  barring  further  testing  and  discovery,  and  survivabil¬ 
ity  assessment  could  proceed  with  these  environments  as  a 
common  input. 


YII.  6 


Apj>enj|lx  Al  lift 
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Overview 

This  appendix  chronicles  the  development  of  NOSVET 
(NOnp a r am e t r i c  'after  SVeeder'  Estimation  Technique),  the 
primary  algorithm  nsed  to  find  the  distribution  of  functions 
of  random  variables.  The  actual  use  of  NOSVET  for  this 
purpose  is  discussed  in  Appendix  C.  In  this  Appendix,  the 
work  of  James  Sweeder  and  A.J.  Moore  [14]  is  examined  as  a 
possible  tool  in  reliability  theory.  One  of  Sweeder's  early 
research  models  (Model  5332)  is  selected  for  further 
development.  The  Model  is  applied  to  random  samples  of  size 
50  from  the  uniform,  normal,  double  exponential,  and  log¬ 
normal  distributions.  These  results  show  unwanted  varia¬ 
tions  due  to  the  randomness  of  the  samples.  A  second  series 
of  experiments  are  performed  to  show  the  performance  of 
Model  5332  on  stylized  samples.  The  concept  of  a  "sty¬ 
lized"  sample  [14]  is  explored  and  defined  more  precisely. 
The  performance  of  the  Model  on  stylized  samples  is 
investigated.  The  results  indicate  a  need  to  change  the 
part  of  the  algorithm  that  affects  the  tails  of  the  distri¬ 
butions.  A  third  series  of  numerical  experiments  is  pre¬ 
sented  that  shows  marked  improvement  in  the  estimate  of  the 
cumulative  distribution  function  (CDF),  but  a  very  much 
degraded  estimate  of  the  probability  density  function  (PDF). 
This  mystery  leads  to  a  fourth  series  of  experiments.  In 


this  series,  Sweeder's  original  trigonometric  interpolation 
of  the  CDF  is  abandoned,  as  is  his  method  of  estimating  the 
PDF.  A  linear  interpolation  scheme  is  implemented,  along 
with  a  centered-dif f erence  scheme  for  finding  the  PDF.  The 
results  are  a  much  improved  estimate  of  both  the  CDF  and  the 
PDF  for  the  distributions  studied.  Finally,  a  fifth  series 
of  experiments  is  performed  in  order  to  find  an  adaptive 
endpoint  extrapolation  technique.  Turning  again  to  numeri¬ 
cal  analysis,  one  finds  that  enforcing  the  conservation  of 
probability  in  the  tails  of  the  distribution  leads  to  opti¬ 
mum  selection  of  the  endpoints.  Sweeder's  extrapolation 
rule  is  shown  to  be  a  special  case  of  the  general  endpoint 
selection  algorithm.  These  results  are  incorporated  in  the 
computer  code  NOSWET,  and  used  in  nuclear  survivability 
assessments . 

Sweeder's  Model  5332 

In  a  recently  published  work  [14],  Sweeder  has 
demonstrated  a  nonp a r am e t r i c  technique  for  estimating 
distribution  and  density  functions.  Sweeder  showed  that, 
given  a  set  of  observations  denoted  by 

{zj},  i=l,m,  where 

then  the  sample  distribution  function  F§(z)  is  defined  by: 

Fs(z)**0  v  z<Zq  (A.l) 

For  the  region  *iiz<zi+1  Fg(z)  is  given  by 

Fs(z)®Gi+[(Gi  +  1-Gi)/2]*[l-cos{n(z-zi)/(zi  +  1-zi)}](A.2) 


Fs(z)*l  v  +  l  (A. 3) 

Sf coder  shoved  that  this  saaiple  distribution  function 
converges  uniforsily  to  the  underlying  distribution  function 

FZU). 

In  Equation  A.2  G ^  is  a  nonparaaetric  plotting  position 
given  by  some  rule  such  as: 

Gj-Ci+aJ/tm+p);  -lloiPil  (A.4) 

In  Equations  A.l  and  A. 3  Zq  and  zm  +  ^  aze  extrapolated 
endpoints  such  that: 

Gq  =  0  (A.  5 ) 

Gm+i=l  (A. 6) 

At  the  data  points  z = z ^  the  sample  distribution 
function  yields  the  values: 

Fs(zi)=Gi  (A. 7) 

Sweeder's  work  vas  examined,  since  it  vas  anticipated 
that  the  problems  of  nuclear  survivability  would  be 
dominated  by  small  sample  statistics,  and  nonparaaetric 
estimation  would  completely  eliminate  the  problems  of 
density  function  identification  and  parameter  estimation. 

Although  the  original  Reference  [14]  should  be  consulted 
for  detail,  Sweeder's  basic  idea  is  to  take  the  sample 
defined  by  the  set  described  above,  and  break  it  into  a 


number  of  subsamples.  Each  subsample  is  treated  as 
representative  of  the  entire  population,  and  Equations  A.l 
through  A. 3  are  used  to  fora  estimates  of  the  CDF.  These 
estimates  are  then  averaged.  A  key  parameter  in  one  of 


Sweeder's  models  is  thus  the  number  of  subsamples  to  take. 
The  first  digit,  5,  from  Model  5332  denotes  the  number  of 
sub  sampl e  s . 

Another  parameter  of  importance  is  the  choice  of  the 
plotting  rule.  Table  A.l  shows  several  selections  that 
follow  under  the  general  description  of  Equation  A. 4.  The 
third  choice,  Hazen's  Rank,  is  the  second  parameter  in  Model 
5332. 

The  third  major  parameter  of  importance  in  Sweeder's 
algorithm  is  the  choice  of  the  extrapolation  constant.  For 
many  of  the  plotting  positions  of  Table  A.l,  the  distribu¬ 
tion  is  not  determined  at  the  endpoints.  That  is,  the  G^'s 
do  not  span  the  entire  space  from  0  to  1.  Thus,  the  sample 
CDF  remains  undefined  in  the  tails  unless  one  chooses  the 


points  z q  and  zm  +  1  in  such  a  way  that 

Fs(z0)=0  (A. 8) 

FS<zm+l)=1  (A’9) 

For  each  subsample,  Sweeder  proposed  a  general  extrapolation 

given  by: 

z0=z1_A1(z2"z1)  (A. 10) 

zm+l=zm+Au(zm_zm-l)  (A. 11) 


where  A^  and  Au  are  the  lower  and  upper  extrapolation  con¬ 
stants.  Table  A. 2  gives  the  choices  of  A^  and  An  that 
Sweeder  examined.  For  Model  5332,  choice  3  was  used. 
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The  CDF,  including  the  endpoints,  can  non  he  con¬ 
structed  using  Sweeder's  technique.  However,  Sweeder  added 
an  additional  feature  in  order  to  reduce  sone  of  the  numeri¬ 
cal  noise  in  the  CDF  estimate.  Once  the  CDF  has  been  deter¬ 
mined,  it  may  be  easily  inverted  by  using  a  Newton-Saphson 
or  other  technique.  Rather  than  inverting  the  distribution 
randomly,  Sweeder  used  the  median  rank  points  (choice  2  of 
Table  A.l)  to  invert  the  distribution.  These  points  were 
then  processed  again  by  the  algorithm,  resulting  in  a  better 
estimate  of  the  CDF.  For  Model  5332,  choice  2  of  Table  A.l 
is  used. 

Model  5332  was  chosen  as  a  starting  point,  since  for  a 
range  of  distributions  examined,  it  behaved  well  with 
respect  to  modified  Cramer  Von  Mises  (CVM)  distance  measures 
[14].  Sweeder  applied  it  to  the  double  exponential,  normal 
and  uniform  distributions,  representing  distributions  with  a 
heavy  tail,  a  moderate  tail,  and  a  short  tail,  respectively. 
Model  5332  was  not  optimum  for  any  of  the  distributions 
considered.  However,  compared  to  the  3  other  models  that 
Sweeder  considered,  it  ranked  second  overall  for  smallest 
CVM  distance  averaged  over  all  three  distributions.  More 
significantly.  Model  5332  was  not  the  worst  choice  for  any 
distribution  considered.  Hence,  it  was  selected  as  the 
starting  point  for  possible  use  in  a  reliability  theory 
approach  to  nuclear  survivability. 


In  the  first  series  of  numerical  experiments,  desig¬ 
nated  as  Series  1,  Model  5332  (vith  two  inversions)  was 
applied  to  the  uniform,  normal,  double  exponential,  and 
lognormal  distributions.  The  lognormal  was  added  as  a  bench 
mark  since  it  is  asymmetric  and  is  popular  vith  other 
survivability  analysts  [49].  Random  samples  of  size  50  were 
drawn  from  each  of  these  distributions,  and  the  CDF  and  PDF 
estimated  by  Sveeder's  method.  The  results  are  illustrated 
in  Figures  A.l  through  A.4.  In  these  figures  the  PDF'S  have 
been  scaled  to  their  peak  values  so  as  to  be  able  to  overlay 
the  PDF  and  CDF  on  the  same  plot.  The  results  for  each 
distribution  are  discussed  briefly  below. 

The  uniform  distribution  shows  a  rather  heavy  tailed 
result  compared  to  the  true  one.  The  true  endpoints  are  at 
0  and  1,  whereas  the  numerical  technique  is  putting  them  at 
-.20  and  1.15.  In  addition,  the  PDF  has  some  oscillations. 
These  oscillations  exceed  the  true  peak  value  of  the  PDF  by 
more  than  50%. 

The  approximation  to  the  normal  is,  at  least  visually, 
somewhat  better.  The  true  endpoints  are  of  course  at  +® 
which  cannot  be  matched  by  this  numerical  algorithm.  A  more 
meaningful  comparison  is  the  98th  and  2nd  percentiles,  which 
should  be  at  1  and  0  respectively.  The  behavior  in  the 
right  hand  tail  is  better  than  that  in  the  left.  The  PDF 
estimate  does  fairly  well  except  for  the  peak  which  should 
be  at  .5.  Here  noise  seems  to  l>e  a  problem. 
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Figure  A. 3.  Model  5332  Applied  To  A  Random  Sample  From 
the  Double  Exponential  Distribution 
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The  results  for  the  double  exponential  (DE)  shown  in 
Figure  A. 3  show  very  clearly  the  concave  shape  of  the  DE 
PDF.  The  peak  of  the  PDF  is  also  Batched  fairly  well,  but 
the  distribution  seeas  to  show  considerable  asymmetry,  being 
shorter  in  the  right-hand  tail  than  in  the  left. 

The  results  for  the  lognormal  are  displayed  in  Figure 
A. 4.  Sweeder's  technique  does  remarkably  well  here. 
Although  not  developed  explicitly  for  asymmetric  density 
functions.  Model  5332  estimates  the  lognormal  CDF  and  PDF 
rather  well.  The  asymmetry  in  the  CDF  is  clearly  evident, 
as  is  the  almost  spiking  PDF. 

Although  the  overall  results  are  favorable,  it  is 
difficult  to  tell  at  this  point  how  useful  Sweeder's  method 
might  be  for  engineering  applications.  In  particular,  it  is 
difficult  to  tell  by  observing  one  random  sample  from  each 
distribution  which  features  are  potential  model  problems  and 
which  are  the  result  of  a  random  draw.  This  question  can  be 
dealt  with  by  considering,  and  defining  in  more  detail,  the 
idea  of  a  stylixed  sample,  and  applying  that  idea  to  a 
second  series  of  numerical  experiments. 

Model  Performance  on  Stylized  Samp  1 e  s 

Sweeder  does  not  strictly  define  the  term  "stylized 
sample'*  [14],  but  basically,  a  stylized  sample  is  one  that 


best  represents  a  given  distribution.  This  can  be  defined 
more  precisely  by  observing  the  following.  Of  all  the 
possible  random  samples  that  one  might  choose,  there  is  one 


particular  sample  that  yields  a  best  estimate  o£  the 
underlying  distribution.  In  this  set,  each  z^  of  the  set 
satisfies: 

Gi“FZ(zi)  (A. 12) 
If  each  ij  so  determined  is  then  plotted  using  the  same 
plotting  rule  by  which  it  was  drawn,  then  by  Equation  A. 7, 
Sweeder's  method  exactly  interpolates  the  distribution 
function  Fz(z)  at  the  data  points.  A  stylized  sample  can 
thus  be  defined  formally: 

Definition:  The  set  of  points  {z*),i=l,m  is  said  to  be 
a  stylized  sample  from  the  population  of  the  random  variable 
Z,  if,  for  every  point  in  the  set, 

FS(zi)=FZ(zi)  (A. 13) 

Stylized  samples  of  size  50  were  drawn  from  each  of  the 
bench  mark  distributions.  The  results  are  displayed  in 
Figures  A. 5  through  A.8  and  are  discussed  below. 

The  results  for  the  uniform  distribution  are  displayed 
in  Figure  A. 5.  Clearly,  the  situation  has  improved.  The 
PDF  shows  very  clear  features  of  the  true  PDF,  and  the  peak 
PDF  value  is  very  close.  The  CDF  is  also  much  improved,  but 
one  still  sees  a  rather  long  tail,  which  is  now  symmetric 
about  .50. 

The  behavior  of  Model  5332  on  a  stylized  sample  from 
the  normal  distribution  is  shown  in  Figure  A. 6.  The  PDF 
estimation  is  almost  outstanding.  However,  one  can  still 
see  that  the  CDF  does  not  match  exactly,  even  though  the 
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The  DE  distribution  as  estimated  by  using  Model  5332  is 
illustrated  in  Figure  A. 7.  Again*  the  PDF  estimate  looks 
very  good,  while  the  CDF  estimate  is  not  quite  as  good. 
Again,  symmetry  about  .SO  is  obvious,  with  slightly  more 
tail  than  in  the  normal  case. 

Finally,  the  lognormal  diatribution  as  estimated  by 
Sweeder's  technique  is  shown  in  Figure  A. 8.  The  PDF 
behavior  is  a  little  hard  to  see  on  such  a  scale  sinoe  it 
rises  so  fast,  but  it  matches  fairly  well.  However,  the  CDF 
seems  to  have  a  noticeable  problem  between  the  70th  and  98th 
percentiles. 

In  summary,  the  Series  2  experiments  allow  one  to  see 
the  power  of  Sweeder's  method  on  stylized  samples.  The 
performance  is  very  good  for  the  normal  and  double 
exponential  distributions,  and  is  less  satisfactory  for  both 
the  uniform  and  lognormal.  The  outstanding  problems  seem  to 
be  undue  weight  in  the  tails  of  the  uniform,  and  somewhat 
weak  estimation  of  the  asymmetric  lognormal.  These 
difficulties  are  overcome  by  actually  modifying  the  model. 

El iminat ion  o  f  Sub  s  amp  1 ing 

At  the  moment,  two  problems  remain:  (a)  the  improper 
tail  weight  given  to  the  uniform  distribution,  and  (b) 
somewhat  poor  performance  in  CDF  estimation,  especially  for 
the  uniform  and  lognormal  distribution  functions. 


These  two  probleas  were  found  to  be  related.  Based  on 


the  previous  definition  of  a  stylized  sample,  if  one  treated 
the  satits  sample  as  a  single  subsample  and  plotted  it  by 
Hazen's  Rule,  then  the  CDF  estimate  should  be  exact  (for 
n*J0)  at  the  CDF  values  .01(.02).99.  The  Series  2  experi¬ 
ments  did  not  yield  exact  results  at  these  points  owing  to 


the  subsampling  and 


tging  process. 


A  third  series  of  experiments  was  performed.  For  this 
series,  the  subsampling  was  eliminated.  Said  another  way, 
Sweeder's  Model  5332  was  altered,  using  his  nomenclature,  to 
Model  1332.  The  results  are  shown  in  Figures  A. 9  through 
A. 12.  For  these  experiments  no  inversions  were  necessary, 
decreasing  computer  execution  time  by  a  factor  of  20  or  so. 

The  use  of  this  new  model  on  a  stylized  sample  from  the 
uniform  distribution  is  shown  in  Figure  A.9.  This  Figure  is 
truly  enlightening,  as  will  be  discussed  shortly.  The  CDF 
estimate  is  nearly  perfect,  and  the  tail  weight  problem  has 
all  but  vanished.  In  short,  the  CDF  estimate  is  precisely 
as  predicted.  However,  the  PDF  estimate  is  much  worse. 
Instead  of  having  a  flat  line  at  the  value  1.0,  sinusoidal 
variations  exist  which  oscillate  between  0  and  n/2. 

The  new  estimate  of  the  normal  distribution  is  dis¬ 
played  in  Figure  A. 10.  Again,  the  CDF  estimate  is  very 


good,  while  the  PDF  estimate  has  failed. 
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The  behavior  of  the  new  model  on  the  DE  and  lognormal 
distributions  is  shown  in  Figures  A. 11  and  A. 12.  Again,  the 
general  result  is  a  very  good  estimate  of  the  CDF  and  a 
badly  oscillating  PDF.  This  problem  led  to  a  fourth  series 
of  numerical  experiments. 

Modification  o f  the  PDF  Es t imat ion  Method 

At  the  conclusion  of  the  Series  3  experiments  two 
observations  were  made.  First,  the  CDF  estimate  was 

behaving  exactly  as  anticipated  based  on  the  concept  of 
stylized  sampling.  Second,  the  PDF  estimate  had  gone  badly 
awry.  Sweeder's  subsampling  technique  clearly  helps  the  PDF 
estimate  a  good  deal.  Why  is  this  so? 

The  PDF  mystery  is  solved  by  examining  in  detail 
Sweeder's  method  for  finding  the  PDF.  He  simply 

differentiated  Equations  A.l  through  A. 3  resulting  in: 

fs(z)=0  v  z<zQ  (A. 14) 

For  the  region  zii.z<zi  +  ^  the  PDF  is  given  by: 
fg(z)  =  (nf  £  +  ^)  s  in ln( z-z  ^ ) / (z  ^  +  z  ^ ) )  (A. IS) 

fg(z)=0  v  z2zm+1  (A. 16) 

The  term  ^i+i/2  in  Equation  A. 15  is  defined  by: 

fi+l/2=<Gi+l"Gi>/<zi+l-2i>  (A'17> 

The  quantity  is  tlie  classic  centered-dif f erence 

approximation  to  the  derivative.  If  the  true  value  of  the 
derivative  is  a  constant  1.0,  one  sees  from  Equation  A. 15 
that  the  PDF  estimate  is  a  sinusoidal  function  with  a  peak 
amplitude  of  n/2,  exactly  as  observed  in  Figure  A. 9. 
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One  also  sees  from  Equation  A. IS  that  increasing  the 
nnaber  of  points  does  nothing  to  iaprove  the  PDF  estimate. 
The  PDF  always  goes  to  zero  at  the  data  points,  and  oscil¬ 
lates  between  0  and  nl  2.  In  short,  the  problea  is  one  of 
posedness,  as  discussed  by  Lee  110]. 

However,  one  can  calculate  the  average  value  of 
Sweeder's  derivative  over  the  interval  [zi*zi+il* 


-S 


■k 


<f0(z) >«  li+l  f_(z)dz/ (z,.-z.) 


S 


i+1 


(A. 18) 


The  result  is  just: 


<f g(z) >*fi+i/2  (A. 19) 

The  average  value  of  Sweeder's  PDF  is  just  the 
numerical  ce nt e r e d-d i f f e r e nee  value.  It  appears  that 
Sweeder's  subsampling  technique  acts  primarily  as  an 
averaging  algorithm. 

If  these  considerations  are  true,  one  can  now  alter 
Sweeder's  method  again.  A  new  CDF  estimate  can  be  defined 
by : 

Fjj(z)=0  v  z<Zq  (A.20) 

On  the  interval 

Fs(z)=Gi  +  (Gi  +  1-Gi) ( z-z j) / (zi  +  1-zi)  (A. 21) 

Fs(z)=l  v  z 2Z„  +  1  (A. 22) 

The  PDF  estimate  is  given  by: 

fs(z)-0  v  z<z0  (A. 23) 


fe(z)“fj/l(l"zft)/(ll/1_zft)  V  ZA <Z  <  Z9  / « 


(A.24) 


On  the  interval  [zj.j/ 2*xi+l/2^: 

£s(*)-f i_i/2+<£i+i/2_£i-l/2) (*~*i-l/2) 

1  (*i+l/2"*i-l/2)  (A*25) 

On  the  interval  lzB-i/2»za+i)  £g(*)  ie  given  by: 

f  s(z>“f*-l/2(zn+l"z)/<z»+l"x*-l/2)  (A.26) 


while 


f«(z)«0  v  z^z. 


(A. 27) 


In  Equation  A. 25  i  ranges  froa  2  to  a-1  and  *i  +  i/2  is 


defined  by: 


*i±l/2“(*i+*i±l)/2 


(A.28) 


The  resnlts  of  the  above  model  on  the  stylized  samples 
froa  the  benchmark  distributions  are  displayed  in  Figures 
A. 13  through  A. 16.  One  final  improvement  will  be 


extrapolation. 


Endpoint 


Even  though  the  results  displayed  in  Figures  A. 13 
through  A. 16  are  very  rewarding  one  nagging  issue  remains — 
that  of  endpoint  extrapolation.  That  is,  how  should  one 
choose  the  points  zQ  and  za+i?  It  is  not  clear  that 
Sweeder's  choice  of  a  constant  and  Au  remain  optimum  for 
the  newly  developed  algorithm.  For  example,  examination  of 
Figures  A. 5  through  A. 8  indicates  that  the  resultant  end¬ 
point  values  do  not  agree  with  those  of  Figures  A. 13  through 
A. 16.  The  uniform  tail  length  is  shorter  in  Figure  A. 13, 
and  the  lognormal  is  longer  on  the  right  and  positive  on  the 


.  •  -V  •. 


A. 27 
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to  a  Double  Exponentially  Distributed  Variable 


left  in  Figure  A. 16,  as  it  should  be.  However,  the  noraal 
and  DB  tails  in  Figures  A. 14  and  A. IS,  r e ape c t ire  1 y,  are 
shorter,  which  they  should  mot  be. 

Sweeder  [14]  grappled  with  this  problem,  developing  a 
number  of  adaptive  models.  This  approach  was  not  completely 
successful,  since  under  subsampling  and  inversion,  endpoints 
were  sometimes  chosen  that  eliminated  some  original  data. 
In  addition,  one  would  like  to  avoid  computer  intensiveness 
as  much  as  possible. 

An  alternate  approach  is  to  enforce  conservation  of 
probability  in  the  tails  of  the  sample  distribution  to  find 
the  endpoints. 

Derivation  of  the  General  Form.  The  numerical  situa¬ 
tion  in  the  tails  is  illustrated  in  Figure  A. 17.  The  CDF  is 
known  at  the  circled  locations.  The  PDF  is  known  approxi¬ 
mately  at  the  location  of  the  x's.  The  locations  of  zQ  and 
zm+^  are  desired. 

The  conservation  of  probability  may  be  approximately 
enforced  on  the  intervals  [zq.z^I  and  [z^,Z2l  by  requiring: 

U-2>) 

<A'3<» 


Using  the  relations: 


j;r  fz<a)dza<fi+i+fi)(*i+i-*i)/2  (A-32> 

Equations  A. 29  and  A.30  rednoe  to: 

Gl“fl(zl"z0)/2  (A.33 ) 

AG-(f2+f1)(z2-z1)/2  (A. 34) 

where 

AGm6i  +  l-6i“G2"Gl  (A.35) 

Since  the  Gj's  are  known  fro*  the  plotting  rule  of 
Equation  A. 4  and  f2  by  some  aethod  of  interpolation  on  the 
interval  tz3/2,z5/2^'  the  only  unknowns  above  are  f2  and  zq. 
Solving  for  Zq  leads  to  the  equation: 

z0-z1-2G1(z2-z1) / {2AG-f 2(z2-z2)  }  (A.3  6) 

If  Aj  is  defined  by: 

A1-2G1/{2AG-f2(z2-z1))  (A.37) 

then  Equation  A.36  has  exactly  the  sane  fora  as  Sweeder's 
extrapolation  rule  Equation  A. 10.  A  similar  analysis  for 
the  right  hand  tail  leads  to  a  similar  form  and  an  extrapo¬ 
lation  constant  Au  defined  by: 

A=2(l-G)  /  {2AG-f  .  (  z  —  z  -  )  }  (A.3  8) 

u  m  a- 1  m  m— l 

Exploiting  the  general  form  for  the  G^'s,  the  extrapo¬ 
lation  constants  may  be  written  as: 

A1-2(l  +  o)/{2-f2(z2-z1)/AG)  (A. 3  9) 

Au-2  (p-a)  /  U-f*.!  (  zm~za_1 )  /AG)  (A.40) 

These  extrapolation  constants  depend  on  the  plotting 
rule  the  sample  size  (m),  the  extreme  order  statis¬ 

tics  (zi*z2'zm-i'zB)*  snd  on  the  choice  of  interpolation  for 


the  derivative  (f~,f _ . ) . 


FORMULA 


2  ( 1+at) 


The  dependence  of  the  eztnpolition  constants  on  the 
plotting  rale  is  best  seen  by  examining  Table  A. 3.  The 
adaptive  extrapolation  constants  properly  reflect  the  under¬ 
lying  behavior  of  the  plotting  rale.  For  exaaple.  the 
eapirical  distribution  function  (EDF)  requires  no  upper 
extrapolation  (Au*0),  while  the  node  rank  requires  no  upper 
or  lower  extrapolation  (A2«0, A^-0). 

One  also  sees  that  Aj  and  An  remain  positive  quantities 
provided : 

2-f2(z2-z1)/A6>0  (A.41) 

2-£m-l(zm-z»-l>M6>0  <A-42> 

Equations  A.41  and  A.42  are  just  the  requirements  that 

f2  and  fm  be  non-negative.  Since  a  numerical  approximation 
is  being  used  for  the  f^'s,  the  above  equations  may  not 

always  be  satisfied.  A  number  of  interpolation  schemes  are 

investigated  below. 

Linear  I  n  t.  e,ip.  °  Ai  t  ion .  From  Figure  A. 17,  and  assuming 

that : 

fO=fm+l“°  (A-43> 

then 

£3/2-^£2-^£5/2  (A. 44) 

fm-l/2^fm-l-^fm-3/2  (A. 45) 

where,  in  general, 

£  i+1  /  2*^®  ^  z  i  +  l~  *  i  ^  (A.46) 

Linear  interpolation  consequently  leads  to: 

f2“f3/2  +  (f5/2"f3/2)  <z2"z3/2)/(z5/2“I3/2) 


(A. 47) 


p. 


I 


i 
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fm-l“fBi-l/2+(fm-3/2"f*-l/2)  (ib-1"z«-1/2) 

1 (zb-3/2_2b-1/2)  (A*48) 
Average  Ve lue  Approx ima  t ion.  The  above  expressions  nay 

lead  to  negative  A^  and  Ag.  Another  choice  for  fj  *a<*  f B_j 


would  be  to  set: 

f  2  =  *  f  3  /  2  + f  5  /  2^2  (A. 49) 

f.-lS(f.-l/2+£.-3/2>/J 
This  choice  leads  to: 

f  2  *  Z2-Z1^A®“*  z3~z1^  C2(*3-Z2>  )  (A. 51) 

*m-l^zm~zm-l^AGm^zm-zm-2^  ^ ^ z m-l_za»-2^  (A. 52) 


In  this  case,  the  resulting  expressions  for  Aj  and  AQ 
remain  positive  provided  that: 

(zj-Zj) / (zj-Xj) <4  (A. 53) 

(zm~zm-2)/(zm-l"zm-2)<4  (A. 54) 

Extreme  Value  Approxima t ion.  The  average  value  approx¬ 
imation  may  still  yield  negative  extrapolation  constants  for 
heavy-tailed  distributions.  Using  an  extreme  value  approxi¬ 


mation,  one  sets: 

*2=*3/2  (A. 55) 

fm-lSfm-l/2  (A. 56) 

The  result  here  is  that: 

f2(z2_zl)M6ml  (A.57) 

fm-l(zm~zm-l) /AG=1  (A. 58) 


This  choice  therefore  guarantees  that  the  extrapolation 
constants  remain  positive.  In  fact,  for  this  approximation: 


Ai “2 ( 1+a) 


(A. 59) 


(A. 60) 


Au=2(p-a) 

For  plotting  under  Hazen's  Rule  (<*=-.5 , 0=0) ,  the  extrapola¬ 
tion  constants  rednce  to  Sveeder's  values  (A^sl,Aa=l). 

Final  Re  snl t  s.  The  final  results  obtained  by  incorpo¬ 
rating  these  new  ideas  are  illustrated  in  Figures  A. 18 
through  A. 22.  Some  discussion  here  is  in  order. 

The  results  for  the  uniform  distribution  are  illus¬ 
trated  in  Figure  A. 18.  It  is  interesting  to  note  that  no 
changes  have  occurred  in  the  endpoints.  This  must  mean  that 
Sveeder's  rule  is  correct.  That  this  is  true  can  be  seen  by 
examining  Equations  A. 47  and  A. 48.  Using  the  fact  that,  for 
the  uniform: 

AG=zi+^-zi  v  i=l,m  (A. 61) 
Equations  A. 47  and  A. 48  reduce  to  Equations  A. 57  and  A. 58. 
The  linear  interpolation  scheme  reduces  naturally  to 
Sveeder's  rule  under  the  assumption  of  a  uniform  distri'  u- 
t  ion . 

The  result  for  the  normal  distribution  is  displayed  in 
Figure  A. 19.  An  increased  tail  length  is  evident  compared 
to  the  plot  of  Figure  A. 14.  An  expanded  left-hand  tail  is 
shovn  in  Figure  A. 20  to  illustrate  the  differences  in  the 
approximations.  In  this  figure  the  smooth  curve  represents 
the  true  PDF,  the  X's  the  adaptive  extrapolation  technique 
just  developed,  and  the  circles  Sveeder’s  rule. 
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NOSWET  Applied  to  a  Stylized  Sample 


The  finel  result  for  the  DE  is  illustrated  in  Figure 
A. 21.  The  sdeptive  technique  has  extended  the  tails  con¬ 
siderably  . 

The  final  result  for  the  lognoraal  is  illustrated  in 
Figure  A. 22.  The  adaptive  endpoint  extrapolation  here 
brings  the  left-hand  tail  very  close  to  zero,  and  extends 
the  right-hand  tail  by  aore  than  a  factor  of  3. 

At  this  point,  the  developaental  task  was  considered  to 
be  finished.  The  coaputer  code  NOSVET  was  written  based  on 
the  above  considerations.  A  listing  is  given  in  Appendix  B, 
and  its  priaary  use  is  discussed  in  Appendix  C. 

Summary 

To  summarize,  Sweeder's  Method  [14]  has  been  Modified 
in  four  ways.  First,  the  concept  of  a  stylized  saaple  was 
defined.  The  perforaance  of  one  of  Sweeder's  early  research 
models  on  stylized  samples  froa  the  uniform,  normal,  double 
exponential,  and  lognormal  distributions  led  to  the  elimi¬ 
nation  of  Sweeder's  subsampling  method.  Second,  Sweeder's 
trigonometric  method  of  interpolating  the  CDF  was  changed  to 
simple  linear  interpolation.  Third,  the  PDF  estimation 
method  was  changed  to  a  centered-dif f erence  scheme.  Fourth, 
a  general  approach  to  endpoint  extrapolation  was  taken  by 
enforcing  the  conservation  of  probability  in  the  tails  of 
the  distribution.  Sweeder's  extrapolation  rule  [14]  was 
seen  to  be  a  special  case  of  this  aore  general  approach. 
The  new  algorithm,  written  as  computer  program  NOSTBT  (NOn- 


v / v.’y.v. /. V. »v* r  •.* k~ r  ■.  ■».  -.  -.  t-.-t 


paraaetric  'after  SVeeder'  Estiaatioa  Technique),  was  found 
to  be  especially  useful  in  finding  the  diatribntion  of 
fsictiois  of  randoa  variables.  The  prograa  listing  is  shown 
in  Appendix  B,  and  its  priaary  use  is  described  in  Appendix 


Appendix  fij.  Listing  Q£  Program  HQSHBT 


10  • ************program  noswet********************** 

20  1  Program  NOSWET  (NOnparametric  'after  SWeeder' 

Estimation  Technique) 

Constructs  a  Nonparametric  Distribution  and  Density 
Function  From  The  Data  in  the  ID  Array  XDATA 

30  . . ............... . 

40  '  Algorithm  is  based  on  NONPARAMETRIC  ESTIMATION  OF 

DISTRIBUTION  AND  DENSITY  FUNCTIONS  WITH  APPLICATIONS 
by  James  Sweeder,  Ph.D,  Capt,  USAF.  AFIT/DS/MA/82-1 

50  . . """" 

60  1  This  Microsoft  Basic  Code  Written  By  HALVOR  A.  UNDEM, 
Capt,  USAF,  DS-83,  as  an  Applications  Tool  in 
Nuclear  Survivability. 

70  ....... ............. ...... 

g0  ************************************************** 

90  'In  the  Fictitious  Fortran  CALL  Statements,  Variables 

Preceding  the  Semicolon  Are  Input,  Those  After  Are 

Output  or  Altered. 

100  .......................... 

110  COMMON  BASICNAME$, SETMIN% , SETMAX% 

120  DIM  MAX(15) ,MIN(15) ,WEIGHT(15) 

130  ON  ERROR  GOTO  5400'  DISK  I/O  ERROR  TRAP 

140  INPUT  "BASICNAME?,SETMIN%,SETMAX%:",BASICNAME?, 

SETMIN% , SETMAX% 

150  INPUT  "Do  You  want  a  DEBUG $  Run  (Y/N)  ", 'DEBUG? 

160  INPUT  " (S)weeder , (H) istogram,  or  (C)ontinuous  PDF"; 
PDFTOGGLE? 

170  PRINT  -THIS  IS  PROGRAM  NOSWET. 015" 

180  '  NOSWET. 015  FEATURES:  (1)— -CHOICE  OF  KSUBS% 

(2) — CHOICE  OF  EXTRAPOLATION 

(3) -- -CHOICE  OF  PDF 

(4)  --WRITES  TO  DISK  IF 

MEMORY  SHORT 

190  '  (5) — CDF  IS  LINEAR 

INTERPOLATION 

200  PREFIX?® "A: " 

210  FOR  SET%=SETMIN%  TO  SETMAX% 

220  PRINT 

230  PRINT  "NOSWET  IS  PROCESSING  DATA  SET  ";SET% 

240  NUMBER?=STR?(SET%) 

250  ADD?sMID? (NUMBER?, 2) 

260  F ILNAM ?-BAS ICNAME $+ ADD ? 

270  'Call  LOADER (; XDATA [] ,SIZE«) 

280  GOSUB  1000 

290  IF  DEBUG?“"Y"  THEN  300  ELSE  370 

300  'THEN  Segment — Debug  Mode  Selected 

310  FOR  I%=1  TO  SIZE% 


320 

330 

340 

350 

360 

370 

380 

390 

400 

410 

420 

430 

440 

450 

460 

470 

480 

490 

500 

510 

520 

530 

540 

550 

560 

570 

580 

590 

600 

610 

620 

630 

640 

650 

660 

670 

680 

690 

700 

710 

720 

730 

740 

750 

760 


PRINT  XDATA(I%); 

IP  1%  MOD  5=0  THEN  PRINT 
NEXT  1% 

PRINT  "SIZEt  IS  "; SIZEt 

IP  SIZE%<2  THEN  PRINT  "SAMPLE  TOO  SMALL- 
LOADER  ABORT": STOP 
•Call  SORT ( XDATA [ ] , N% ; XDATA [ ] ) 

N%=SIZE% 

GOSUB  1490 

•Call  SUBSAMPL ( XDATA [ ] ,KSUBS% , SIZE% , ; SUBSAMP [ ] ,M%,R%, 
XMAX, XMIN) 

INPUT  "(R)andom  or  (S)tylized  Sample"; TOGGLE $ 

IP  TOGGLE$-"R"  THEN  KSUBSt=5  ELSE  KSUBS%=1 
IF  KSUBS%  >N%/3  THEN  PRINT  "RSUBSt  TOO  LARGE" :GOTO 
410 

INPUT  "(S)weeder  or  (A) utoranging  Extrapolation"; 

EXTRAPTOGGLE  $ 

GOSUB  1780 

IF  DEBUG$«"Y*  THEN  470  ELSE  580 
•THEN  Segment — Debug  Mode  Selected 
FOR  SAMPLEt =1  TO  RSUBSt 
IP  SAMPLEt <=Rt  MEN  LASTELfMt  +  2  ELSE 
LASTELt  *M%  4-1 

PRINT  "SAMPLEt ,  LASTELt  ARE  ";  SAMPLEt ;  LASTELt 

FOR  ELEMNTt-0  TO  LASTELt 

PRINT  SUBSAMP ( ELEMNTt , SAMPLEt ) ; 

IF  ELEMNTt  MOD  5=0  THEN  PRINT 
NEXT  ELEMNTt 

PRINT  "PAUSING  BEFORE  NEXT  SAMPLE": STOP 
NEXT  SAMPLEt 

PRINT  "Mt,Rt, XMAX, XMIN  ARE  ”;Mt;Rt; XMAX; XMIN 
•Call  JACKNIFE( SUBSAMP l ] , SIZEt , RSUBSt ; SUBSAMP [ ] f 
XDATA [] , ZDATA [ ] ) 

GOSUB  2600 

IF  DEBUG$="Y"  THEN  610  ELSE  820 
•THEN  Segment — Debug  Mode  Selected 

PRINT  "HERE  IS  THE  NEW  XDATA  ARRAY"; 

FOR  I%=1  TO  SIZEt 
PRINT  XDATA (It); 

IF  1%  MOD  5=0  THEN  PRINT 
NEXT  1% 

PRINT  "HERE  IS  THE  ZDATA  ARRAY" 

FOR  I%=1  TO  SIZEt 
PRINT  ZDATA (It); 

IF  It  MOD  5=0  THEN  PRINT 
NEXT  It 

PRINT  "HERE  IS  THE  NEW  SUBSAMP  ARRAY" 

FOR  SAMPLEt =1  TO  RSUBSt 
IF  SAMPLEt <=Rt  THEN  LASTELt =Mt+2  ELSE 
LASTELt =Mt+l 

PRINT  "SAMPLEt, LASTELt  ARE  "; SAMPLEt ; LASTELt 
FOR  ELEMNTt =0  TO  LASTELt 


770  PRINT  SUBSAMP ( ELEMNT% , SAMPLE% ) ; 

780  IF  ELEMNT%  MOD  5=0  THEN  PRINT 

790  NEXT  ELEMNT% 

800  PRINT  "PAUSING  AFTER  SAMPLE%  " ; SAMPLE% : STOP 

810  NEXT  SAMPLE% 

820  INPUT  "Dump  XDATA  To  Disk  ";ANS$ 

830  IF  ANS$="Y"  THEN  840  ELSE  860 
840  'THEN  Segment — Dump  XDATA  To  Disk 
850  GOSUB  5160 

860  INPUT  "Want  A  Look  (Y/N)";LOOK$ 

870  IF  LOOK$="Y"  THEN  880  ELSE  940 
880  INPUT  "Input  The  Value  of  X  ",X 

890  'Call  CDFPDF( SUBSAMP [] ,KSUBS%,M%,R%, IPLOT% ,X; AVGCDF (X) , 

AVGPDF (X) ) 

900  IPLOT%=3  'Midpoint  of  EDF 

910  GOSUB  3780 

920  PRINT  "AT  X=";X; "CDF, PDF  ARE  "; AVGCDF ; AVGPDF 
930  GOTO  860 

940  INPUT  "Do  You  Wish  To  Plot  the  PDF/CDF  ";ANS$ 

950  IF  ANS$="Y"  THEN  GOSUB  5450'  CALL  PLOTPCDF 

960  NEXT  SET% 

970  CHAIN  "FAILPROB" 

980  END 

990  '""""""" . 

1000  SUBROUTINE  LOADER ( ,•  XDATA [], SIZE%)  """"  ' 

1010 

1020  'Subroutine  LOADER  Loads  the  XDATA  Array  From 
The  Keyboard  or  From  A  Datafile 

1030  . """' 

1040  'input  VARIABLES:  None— Prompts  For  All 
1050 

1060  'OUTPUT  VARIABLES:  XDATA — The  Array  Containing  SIZE% 

Elements 

1070 

1080  IF  DEBUG$="Y"  THEN  PRINT  "LOADER  HAS  BEEN  CALLED" 

1090  ANS$="N" 

1100  IF  ANS$="Y"  THEN  1110  ELSE  1190 

1110  'THEN  Segment — Load  XDATA  From  Random  Number  Generator 
1120  RANDOMIZE 

1130  INPUT  "What  Size  is  Your  Sample  ";SIZE% 

1140  DIM  XDATA (SIZE%) 

1150  FOR  J%=1  TO  SIZE% 

1160  XDATA (J%)=RND 

1170  NEXT  J% 

1180  RETURN 

1190  'ELSE  Segment — Load  From  Keyboard  or  Datafile 
1200  KEYIN$="T" 

1210  IF  KEYIN$="K"  THEN  1220  ELSE  1290 

1220  'THEN  Segment — Load  Data  From  Keyboard 

1230  INPUT  "What  Size  Is  Your  Sample" ;SIZE% 

1240  DIM  XDATA(SIZE%) 

1250  FOR  J%=1  TO  SIZE% 
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1260 

1270 

1280 

1290 

1300 

1310 

1320 

1330 

1340 

1350 

1360 

1370 

1380 

1390 

1400 

1410 

1420 

1430 

1440 

1450 

1460 

1470 

1480 

1490 

1500 


PRINT  "XDATA(";J!; ")=";: INPUT  XDATA(J%) 
NEXT  J% 

RETURN 

1  ELSE  Segment — Load  Data  From  Tape 
OPEN  ■I",3,FILNAM$ 

INPUT  #3  f  SIZE% t NEL% 

IF  TAPCALL$="T"  THEN  ERASE  XDATA 
DIM  XDATA <SIZE%) 

NREC%  =  S I Z  E% \NEL% 

IF  SIZE%  MOD  NEL%<>0  THEN  NREC!=NREC!+1 
XINDEX%=0 

FOR  RECORD! =1  TO  NREC% 

FOR  ELEMNT%=1  TO  NEL« 

X INDEX! =XINDEX%+1 

IF  XINDEX!>SIZE!  THEN  PRINT  "OUT  OF 
DATA": RETURN 
INPUT  #3 , XDATA (XINDEX!) 

NEXT  ELEMNT! 

NEXT  RECORD! 

CLOSE  *3 
TAPC ALL  $ = " T " 

RETURN 

END 


'  SUBROUTINE  SORT ( XDATA [ ] , N! ; XDATA [])«■"■» 

I nnnnnnnnnnnnnnnnnnnnnnnnnnn 


1510  'SUBROUTINE  SORT  Sorts  a  Data  Set  From  Min  to  Max 
1520  'INPUT  VARIABLES:  XDATA — The  Array  Containing  the  Data 

N% - The  Number  of  Elements 


1530  'OUTPUT  VARIABLES:  XDATA— The  Array  After  Sorting 

1540  ' ********************************************** 

1550  IF  DEBUG$="Y"  THEN  PRINT  "SORT  HAS  BEEN  CALLED" 

1560  FLIPS=1  'FORCE  AT  LEAST  ONE  PASS 

1570  WHILE  FLIPS 

1580  FLIPS=0 

1590  FOR  J!=l  TO  N!-l 

1600  IF  XDATA (J! ) >XDATA(J!+1)  THEN  1620  ELSE  1660 

1610  'THEN  Segment — SWAP  pair 

1620  SWAP  XDATA(J!) fXDATA(J!+l) 

1630  FLIPS® 1 

1640  GOTO  1660 

1650  'ELSE  Segment — Look  At  Next  Pair 

1660  NEXT  J! 

1670  WEND 

1680  IF  DEBUGS® "Y"  THEN  1690  ELSE  1750 
1690  PRINT  "HERE  IS  THE  SORTED  ARRAY" 

1700  FOR  I%=1  TO  N! 

1710  PRINT  XDATA (I!) ; 

1720  IF  I!  MOD  5*0  THEN  PRINT 
1730  NEXT  I! 

1740  PRINT 
1750  RETURN 
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1760  END 

1770  . . ■■"""■ 

1780  1  SUBROUTINE  SUBSAMPL ( XDATA [ ] , KSUBS% f SIZE% ; SUBSAMP [ ] , 

M% , R% f XMAXr XMIN) 

1790  . . . 

1800  'Subroutine  SUBSAMPL  Loads  the  2D  Array  SUBSAMP  from 
the  ID  Array  XDATA 

1810 

1820  'INPUT  VARIABLES;  XDATA — ID  Array  Containing 

Observations 

SIZE% — Total  Dimension  of  XDATA 
KSUBS%-Number  of  Subsamples 
Desired — Will  Be  #  of  Cols  of  SUBSAMP 
1830  !■«•■■■■■■■■■■■■■•■■■■■■■■■ 

1840  'OUTPUT  VARIABLES:  M% — Nominal  Number  of  Elements  Per 

Subsample 

R% — Number  of  Subsamples  With  M%+1 
Elements 

SUBSAMP — 2D  Array — Each  Column  Is  A 
Subsample 

1850 

1860  'REQUIRED  EXTERNALS:  SUBROUTINE  ENDPOINT — To  Get 

XMIN,XMAX 

1870 

1880  IP  DEBUG$=*Y"  THEN  PRINT  "SUBSAMPL  HAS  BEEN  CALLED" 

1890  M%=SIZE%\KSUBS% 

1900  R%=SIZE%  MOD  KSUBS% 

1910  IF  S AMPC ALL% * 1  THEN  1920  ELSE  1940 
1920  'THEN  Segment— -SUBSAMPL  Previously  Called 
1930  ERASE  SUBSAMP 

1940  IP  R%=0  THEN  1950  ELSE  1980 

1950  'THEN  Segment — M%  Elements  in  SUBSAMP — Dimension  1  More 
1960  DIM  SUBSAMP (M%+1,KSUBS%) 

1970  GOTO  2000 

1980  'ELSE  Segment — M%+1  Elements  in  SUBSAMP 
1990  DIM  SUBSAMP (M%+2 , KSUBS%) 

2000  XINDEX%=0 

2010  FOR  SAMPLE%=1  TO  KSUBS% 

2020  IF  SAMPLE%<=R%  THEN  2030  ELSE  2060 

2030  'THEN  Segment — M%+1  In  This  One 

2040  SAMPSIZE%=M%+1 

2050  GOTO  2080 

2060  'ELSE  Segment — M%  In  This  One 

2070  SAMPSIZE%=M% 

2080  WEIGHT (SAMPLE%) =SAMPSIZE%/SIZE% 

2090  FOR  ELEMNT%=1  TO  SAMPSIZE% 

2100  XINDEX%=XINDEX%+1 

2110  IF  TOGGLE $="R"  THEN  XINDEX%=SAMPLE%+KSUBS%* (ELEMNT%-1) 
2120  SUBSAMP ( ELEMNT% , SAMPLE% ) =XDATA ( X INDEX% ) 

2130  NEXT  ELEMNT% 

2140  'CALL  EXTRAP ( SUBSAMP [ ] , TOGGLE? j DELTAL , DELTAU ) 

2150  GOSUB  5970 


2160  SUBSAMP ( 0 , SAMPLE% ) =SUBSAMP ( 1 , SAMPLE% ) -DELTAL 
* ( SUBSAMP ( 2 , SAMPLE % ) - SUBSAMP ( 1 , SAMPLE % ) ) 

2170  SUB  SAMP ( SAMPS IZE%+1 , SAMPLE! ) = SUBSAMP ( SAMPS IZE! , SAMPLE! ) 
+DELTAU* ( SUBSAMP ( SAMPS IZ  E! , SAMPLE! ) 

-SUBSAMP (SAMPSIZE!-1 f SAMPLE!) ) 

21 80  MIN ( SAMPLE! ) = SUBSAMP ( 0 , SAMPLE! ) 

2190  MAX ( SAMPLE! ) = SUB SAMP ( SAMPS IZ E! +1 , SAMPLE! ) 

2200  NEXT  SAMPLE! 

2210  'CALL  SUBROUTINE  ENDPOINT (MIN [] ,MAX[] ,KSUBS!) 

2220  GOSUB  2270 
2230  SAMPCALL! =1 
2240  RETURN 
2250  END 

2260  . . nmnnmnmumn, 

2270  '  SUBROUTINE  ENDPOINT (MIN [] , MAX [] , KSUBS! ?XMIN,XMAX) 
2280 

2290  'Subroutine  ENDPOINTS  Gets  The  Extrapolated  Values 
for  A  Set  of  Data  Found  in  Arrays  MIN  And  MAX 

2300 

2310  'INPUT  VARIABLES:  MIN — ID  Array  Containing  Mininura 

Ex traps 

MAX — ID  Array  Containing  Max  Extraps 
RSUBS! — Size  of  Above  Arrays 

2320 

2330  'OUTPUT  VARIABLES:  XMIN — Minimum  Found  in  MIN 

XMAX — Maximum  Found  in  MAX 

2340  .  iM.miH...... 

2350  'REQUIRED  EXTERNALS:  SUBROUTINE  SORT 

2360  . . """""""' 

2370  IF  DEBUG $="Y"  THEN  PRINT  "ENDPOINT  HAS  BEEN  CALLED" 

2380  ' CALL  SORT ( XDATA [ ] , N! ; XDATA [ ] ) 

2390  FOR  J!=l  TO  KSUBS! 

2400  SWAP  XDATA (J!) f MIN (J!) 

2410  NEXT  J! 

2420  N! -KSUBS! 

2430  GOSUB  1490 

2440  FOR  J!=l  TO  KSUBS! 

2450  SWAP  XDATA (J!) rMIN(J!) 

2460  NEXT  J! 

2470  XMIN-MIN(l) 

2480  ' CALL  SORT (XDATA [] , N! ; XDATA [ ] ) 

2490  FOR  J!=l  TO  KSUBS! 

2500  SWAP  XDATA (J!) r MAX (J!) 

2510  NEXT  J! 

2520  GOSUB  1490 

2530  FOR  J!=l  TO  KSUBS! 

2540  SWAP  XDATA (J!) rMAX(J!) 

2550  NEXT  J! 

2560  XMAX=MAX (KSUBS!) 

2570  RETURN 
2580  END 


2600  '  SUBROUTINE  JACKNIFE (SUBSAMP [ ] , SIZE% ,RSUBS% ; 

SUBSAMP I ] , XDATA [ ] , ZDATA [ ] ) 

2610  •""»»»"*  . . 

2620  'Subroutine  JACKNIFE  Samples  From  the  Current 
Distribution,  Storing  the  Sample  in  ZDATA 
2630  »■■■■■■■■■■■■■■■■■■■■■•■■ 

2640  'INPUT  VARIABLES:  SUBSAMP — 2D  Array  of  Subsamples  For 

Current  Distribution 

SIZE% - Total  Number  of  Data  Points 

in  Original  Sample 

KSUBS% — Number  Of  Subsamples  Desired 
2650  . . ""• 

2660  'OUTPUT  VARIABLES:  ZDATA — Array  Containing  Last 

Pseudodata 

SUBSAMP — Altered  2D  Array  For  "New" 
Distribution 

XDATA — Now  Contains  New  Dataset 

2670 

2680  'REQUIRED  EXTERNALS:  SUBROUTINE  PLOTPNT — To  Get  Median 

Ranks 

SUBROUTINE  ZOFCDF — To  Invert 

Distribution 

SUBROUTINE  SUBSAMPL — To  Reload  Array 

SUBSAMP 

2690 

2700  IF  DEBUGS— "Y"  THEN  PRINT  "JACKNIFE  HAS  BEEN  CALLED" 

2710  IF  JACKCALL$<>"T"  THEN  DIM  ZDATA(SIZE%) 

2720  INPUT  "Input  JACKMAX%  (0,1,2) : ",JACKMAX% 

2730  IF  JACKMAX%=0  THEN  JACKCALL$«"T" : RETURN 
2740  IF  JACKCALL$0"T"  AND  JACKMAXOO  AND  T0GGLE$="S"  THEN 
2750  ELSE  2780 

2750  'THEN  SEGMENT — Verify  Desired  Operation 
2760  PRINT  "YOU  ARE  JACKNIFING  A  STYLIZED  SAMPLE" 

2770  STOP 

2780  FOR  JACK%=1  TO  JACKMAX% 

2790  PRINT 

2800  PRINT  "THIS  IS  INVERSION  NUMBER  ";JACK% 

2810  ZDATA (0)=XMIN 

2820  FOR  IJACK%=1  TO  SIZE% 

2830  ' CALL  PLOTPNT { IPLOT% , IPNT% , NOBS% ; PLOTPNT) 

2840  IPNT%=IJACK% 

2850  IPLOT%=2 '  Median  Rank  Points 

2860  NOBS%sSIZE% 

2870  GOSUB  3060 

2880  'CALL  ZOFCDF (ALFA, ZMIN, ZMAX; ZALFA) 

2890  ZMAXsXMAX 

2900  ZMIN= ZDATA (IJACK%-1) 

2910  ALFA= PLOTPNT 

2920  GOSUB  3270 

2930  ZDATA (IJACK%)=Z ALFA 


2940  PRINT  "FOR  DATA  POINT  IJACK%; "ALFA, ZDATA  ARE  "jALFA 
Z  DATA ( I JACK% ) 

2950  NEXT  IJACK% 

2960  1  CALL  SUBSAMPL(XDATA[] , KSUBS% , SIZE% ; SUBSAMPL [ ] ,M%,R%, 

XMAXf XHIN) 

2970  FOR  J%=1  TO  SIZE% 

2980  SWAP  XDATA(JS) , ZDATA (J%) 

2990  NEXT  J% 

3000  GOSUB  1780 

3010  NEXT  JACK% 

3020  JACKCALL$*"T" 

3030  RETURN 
3040  END 
3050 

3060  '  SUBROUTINE  PLOTPNT ( IPLOT% , IPNT% , NOBS% ; PLOTPNT) 

3070  i 

3080  1  Subroutine  PLOTPNT  Generates  the  Plotting  Position 
for  the  I%th  Observation  out  of  NOBS% 

3090  . . . . 

3100  'INPUT  VARIABLES:  IPLOT% — Variable  Selecting  Choice 

of  Plotting  Position 

IPNT% - The  I%th  Observation 

NOBS% - The  Total  Number  of 

Observations 

3110  . . 

3120  'OUTPUT  VARIABLE:  PLOTPNT — The  Plotting  Position 
3130  . """"' 

3140  IF  DEBUG $="Y"  THEN  PRINT  "PLOTPNT  HAS  BEEN  CALLED" 

3150  IF  IPNT%=0  THEN  PLOTPNT* 0 : RETURN 

3160  IF  IPNT%=NOBS%+l  THEN  PLOTPNT*l : RETURN 

3170  IF  IPLOT%*2  THEN  3180  ELSE  3210 

3180  'THEN  Segment — Approximate  Median  Ranks 

3190  PLOTPNT* ( IPNT%- . 3) / (NOBS%+. 4) 

3200  GOTO  3230 

3210  'ELSE  Segment — Midpoint  of  EDF 
3220  PLOTPNT* ( IPNT%~ . 5 ) /NOBS% 

3230  IF  PLOTPNT<0  OR  PLOTPNT>l  THEN  PRINT  "ERROR  IN 
PLOTPNT": STOP 
3240  RETURN 
3250  END 

3260  . . . "*"" 

3270  '  SUBROUTINE  ZOFCDF (ALFA, ZMIN, ZMAX; ZALFA) 

3280  •■■■■■■■■■■■■■■■■■■■■■■■■■■■ 

3290  'Subroutine  ZOFCDF  Finds  the  Value  of  z  such  that 
Pr [ Z<=z ] =ALFA 

3300  •■■■■■■■■■■■■■■■■■■■■■■■■« 

3310  'INPUT  VARIABLES:  ALFA— PERCENTILE  DESIRED 

ZMIN— MINIMUM  VALUE  OF  VARIATE 
ZMAX— MAXIMUM  VALUE  OF  VARIATE 
3320  . . """""""" 

3330  'OUTPUT  VARIABLES:  ZALFA— Value  of  Z  Satisfying 
Equation 


3340 

3350 

3360 

3370 

3380 

3390 

3400 

3410 

3420 

3430 

3440 

3450 

3460 

3470 

3480 

3490 

3500 

3510 

3520 

3530 

3540 

3550 

3560 

3570 

3580 

3590 

3600 

3610 

3620 

3630 

3640 

3650 

3660 

3670 

3680 

3690 

3700 

3710 

3720 

3730 

3740 

3750 

3760 

3770 

3780 


•REQUIRED  EXTERNALS:  SUBROUTINE  CDFPDP— To  Supply 

CDF  Values 


IF  DEBUG $="Y"  THEN  PRINT  "ZOFCDF  HAS  BEEN  CALLED" 

ICOUNT%=0 

CRIT=1 

HIVALUE=ZMAX 
LOVALUE=  ZMIN 
WHILE  CRIT>. 000001 

IF  HIVALUE=LOVALUE  THEN  3440  ELSE  3480 
•THEN  Segment — Nonconvergence  Problem 
PRINT  "ZOFCDF  CANNOT  CONVERGE" 

PRINT  "ALFA, Z ALFA, CRIT  ARE  "; ALFA; ZALFA;CRIT 
STOP 

'ELSE  Segment — Normal  Search  Continues 
ICOUNT% = ICOUNT% +1 

IF  ICOUNT%>212  THEN  3510  ELSE  3580 
"THEN  Segment — Tolerance  Not  Met 
PRINT  "TOLERANCE  NOT  MET" 

PRINT  "ALFA, Z ALFA, CRIT  ARE  "; ALFA; GUESS; CRIT 
PRINT  "HIVALUE,LOVALUE  ARE  " ; HIVALUE; LOVALUE 
INPUT  "Execute  Recovery  Routine  ";ANS$ 

IF  ANS$="Y"  THEN  X=(HIVALUE+LOVALUE)/2:GOTO 
3740 

PRINT  "ZOFCDF  HAS  ABORTED" : STOP 
IF  PDFO0  AND  IC0UNT%>1  THEN  3590  ELSE  3620 
'THEN  Segment — Newton-Raphson  Estimate 
GUESS=GUESS+ (ALFA-CDF) /PDF 
IF  GUESS>HIVALUE  OR  GUESS<LOVALUE  THEN  3630 
ELSE  3640 

'ELSE  Segment — Halve  the  Interval  Estimate 
GUESS- (HIVALUE+LOVALUE) /2 
' CALL  CDFPDF { SUBSAMP  £ ] , KSUBS% , M% , R% , IPLOT% , X ; 

AVGCDF (X) , AVGPDF (X) ) 

IPLOT%=3  'Midpoint  of  EDF 

X=GUESS 

GOSUB  3780 

CDF=AVGCDF 

PDF*AVGPDF 

IF  CDF > ALFA  THEN  HIVALUE=GUESS 
IF  CDF<ALFA  THEN  LOVALUE=GUESS 
CRIT=ABS (CDF- ALFA) 

WEND 

ZALFA=X 

RETURN 

END 


'  SUBROUTINE  CDFPDF ( SUBSAMP I J , KSUBS% , M% , R% , IPLOT% , X ; 

AVGCDF (X) , AVGPDF (X) ) 


>-A 

A 


I 


3790  ' 


I  "  W.J 


3800  'Subroutine  CDFPDF  Gets  the  Value  of  The  CDF  and 
PDF  Using  Sweeder's  Estimation  Technique 

3810 

3820  'INPUT  VARIABLES:  KSUBS%~THE  NUMBER  OF  SUBSAMPLES 

X - VARIATE  OF  DISTRIBUTION 

SUBSAMP-2D  ARRAY  CONTAINING 

PARTITIONED  SAMPLE 
R%— NUMBER  OF  COLS  OF  SUBSAMP  WITH 
M%+1  ELEMENTS 

M% --NUMBER  OF  ELEMNTS  IN  SAMPLES% 
R%+1  ON 

3  830  '  IPLOT%— CHOICE  OF  PLOTTING  POSITION 

3840 

3850  'OUTPUT  VARIABLES:  AVGCDF — Average  CDF  Over  KSUBS% 

Subsamples 

AVGPDF — Average  PDF  Over  KSUBS% 
Subsamples  (Both  Evaluated  at  X) 

3860 

3870  'EXTERNALS  REQUIRED:  SUBROUTINE  LOADXDUM — To  Load 

Dummy  Array 

SUBROUTINE  POINTCDF— To  Get 
Point  CDF , PDF  Valus 

3880  '"""""" . 

3890  IF  DEBUGS* "Y"  THEN  PRINT  "CDFPDF  HAS  BEEN  CALLED" 

3900  SUMCDF=0 
3910  SUMPDF=0 

3920  FOR  SAMPLE%*1  TO  KSUBS% 

3930  'CALL  LOADXDUM (SUBSAMP [J ,XDUM [] , SAMPLE% rR% ,M%) 

3940  GOSUB  4080 

3950  'CALL  POINTCDF (XDUM[ ] , SAMPSIZE% ,X, IPLOT%) 

3960  SAMPS  I Z  E%  *LASTEL% 

3970  GOSUB  4430 

3980  IF  DEBUG$="Y"  THEN  3990  ELSE  4000 

3990  PRINT  "FOR  SAMPLE* " ; SAMPLE% ; "CDF , PDF  ARE  " ; CDF , PDF 
4000  SUMCDF=SUMCDF+WEIGHT ( SAMPLE% ) *CDF 
4010  SUMPDF=SUMPDF+WEIGHT ( SAMPLE! ) *PDF 
4020  NEXT  SAMPLE« 

4030  AVGCDF* SUMCDF 
4040  AVGPDF=SUMPDF 
4050  RETURN 
4060  END 

4070  i  ■■■■■■■■■■■■■■■■■■■■■■■■■■■ 

4080  '  SUBROUTINE  LOADXDUM ( SUBSAMP [  ] ,  M%  ,  R%  , SAMPLE % ; XDUM [ J ) 
4090  . . . 

4100  'Subroutine  LOADXDUM  Loads  the  XDUM  Array  For  Later 
Use  By  Subroutine  POINTCDF — It  Acts  As  The 
Fictitious  Data  Array  And  is  Loaded  From  Array  SUBSAMP 
4110  . """"""' 


B.10 


4120  'INPUT  VARIABLES:  SAMPLE%— IDENTIFIES  COLUMN  OF 

SUBSAMP 

R% - FIRST  R%  COLUMNS  OF  SUBSAMP 

HAVE  ONE  MORE  DATA  Point 
SUBSAMP — 2X2  Array  Containing  Data — 
ROW  is  ELEMNT% ,  Col  is  SAMPLE% 

M% - Number  of  Elements  Per 

Sample 

4130  . . . . 

4140  'OUTPUT  VARIABLES:  XDUM - DUMMY  SORTED  DATA  Array-lD 

4150 

4160  IF  DEBUG$="Y"  THEN  PRINT  "LOADXDUM  HAS  BEEN  CALLED" 

4170  IF  DEBUG?- "Y"  THEN  PRINT  " SAMPLE% , R% , M%  CAME  IN  AS 
" ; SAMPLE% ;R% ,M% 

4180  IF  SAMPLE%<=R%  THEN  4190  ELSE  4220 

4190  'THEN  Segment — Set  Has  M%+1  Elements 

4200  LASTEL% =M% +2  'Because  of  Extrapolated 

Points 

4210  GOTO  4240 

4220  'ELSE  Segment — Set  Has  M%  Elements 
4230  LASTEL% = M%  +1 

4240  IF  CALLXDUM$="T"  THEN  ERASE  XDUM 
4250  DIM  XDUM ( LASTEL% ) 

4260  FOR  I%=0  TO  LASTEL% 

4270  XDUM ( 1% ) -SUBSAMP ( 1%  *  SAMPLE % ) 

4280  NEXT  I« 

4290  IF  DEBUG$="Y"  THEN  4300  ELSE  4360 
4300  'THEN  Segment — Debug  Mode  Selected 

4310  PRINT  "HERE  IS  THE  XDUM  ARRAY  FOR  SAMPLE  SAMPLE* 
4320  FOR  IDEX%=0  TO  LASTEL% 

4330  PRINT  XDUM (IDEX%) ; 

4340  IF  IDEX%  MOD  5=0  THEN  PRINT 

4350  NEXT  IDEX* 

4360  'ELSE  Segment — Normal  Termination 
4370  CALLXDUM$="T" 

4380  RETURN 

4390  END 

4400  . . nwmnn*mnmmmn*wnm.mn 

4410  'SUBROUTINE  POINTCDF (XDUM [ ] , SAMPSIZE* , Xf IPLOT% ; 

CDF, PDF) 

4420  »■■■«•■■■■■■■■■■■■■■■■■■■■■■■ 

4430  '  Subroutine  POINTCDF  Estimates  the  CDF  and  PDF 
Nonparametrically  From  a  Data  Set  Found  in  XDUM — 

4440 

4450  'INPUT  VARIABLES:  XDUM— DUMMY  ARRAY  CONTAINING  THE 

DATA  INCLUDING  THE  ENDPOINTS 
SAMPS  I  ZE%— TOTAL  SIZE  OF  THE  ARRAY 

X - POINT  WHERE  CDF, PDF  WANTED 

4460  '  IPLOT% — Choice  of  Plotting  Position 

4470 

4480  'OUTPUT  VARIABLES:  CDF— Value  Of  CDF  At  Point  X 

PDF— Value  of  PDF  At  Point  X 


SUBROUTINE  PLOTPNT 


4490  .................... 

4500  'REQUIRED  EXTERNALS: 

4510  .................... 

4520  IF  DEBUGS* "Y"  THEN  PRINT  "POINTCDF  HAS  BEEN  CALLED" 
4530  PI-4*ATN(1) 

4540  DUMIN»XDUM(0) 

4550  DUMAX*XDUM ( SAMPS IZ E% ) 

4560  IF  X<-DUMIN  THEN  CDF*0: PDF* 0: RETURN 
4570  IF  X>=DUMAX  THEN  CDF=1 : PDF=0 : RETURN 
4580  FOR  I%*0  TO  SAMPSIZE%-1 

4590  IF  XDUM(I%) <*X  AND  X<XDUM(I%+1)  THEN  4600  ELSE  4710 
4600  'THEN  Segment — Interval  Found 


4610 

'CALL 

PLOTPNT ( 1% , IPLOT% , SAMPSIZE% ) 

4620 

IPNT%® 1% 

4630 

NOBS% = SAMPS I Z  E« -1 

4640 

GOSUB  3060 

4650 

G* PLOTPNT 

4660 

'CALL 

PLOTPNT ( I%+1 , IPLOT% f  SAMPSIZE% 

4670 

IPNT%=I%+1 

4680 

GOSUB  3060 

4690 

GPLUS* PLOTPNT 

4700 

GOTO 

4740 

4710  'ELSE  Segment — Look  Again 
4720  NEXT  1% 

4730  PRINT  "ABORT  IN  POINTCDF — X  NOT  FOUND": STOP 
4740  ARG= (X-XDUM{ 1%) )/{XDUM(I%+l) -XDUM(I%) ) 

4750  CDF=G+ (GPLUS-G) *ARG 

4760  IF  PDFTOGGLE$*"S"  THEN  CDF*G+  (GPLUS-G)  /2 

* (1-COS (PI*ARG) ) 

4770  IF  CDF<-6E-09  OR  CDF>1  THEN  PRINT  "CDF  ABORT  IN 
POINTCDF": STOP 

4780  IF  PDFFLAG$="S"  THEN  PRINT  "PDF  ABORTED " : PDF*0 : RETURN 
4790  IF  PDFTOGGLE$="S"  THEN  4800  EL;  4830 
4800  '  THEN  Segment — Use  Sweeder's  (  iginal  PDF 
4810  PDF=PI/2* (GPLUS-G) / (XDUM ( 1% . 1) -XDUM ( 1% ) ) *SIN (PI*ARG) 
4820  GOTO  5130 

4830  IF  PDFTOGGLE$="H"  THEN  PDF* (GPLUS-G) /(XDUM (I%+1) - 
XDUM(I%) ) :GOTO  5130 

4840  IF  X<=.5* (XDUM(I%) +XDUM(I%+1) )  THEN  4850  ELSE  5020 
4850  'THEN  Segment — Interpolate  From  Previous  Derivative 
4860  IF  XDUM(I%+1)  *XDUM(I%)  THEN  4870  ELSE  4890 
4870  PRINT  "DUPLICATE  SAMPLE— PDF  ABORT" 

4880  PDFFLAG$«"S": PDF-0: RETURN 

4890  PDFMAX* (GPLUS-G) / (XDUM ( I%+1) -XDUM ( 1%)  ) 

4900  XHI=.5*(XDUM(I%)+XDUM(I%+1) ) 

4910  IF  I%*0  THEN  PDFMIN*0 :XLO*DUMIN:GOTO  5120 

4920  'CALL  PLOTPNT ( I%-1 ) fIPLOT%,SAMPSIZE%; PLOTPNT) 

4930  IPNT%=I%-1 

4940  GOSUB  3060 

4950  GMINUS* PLOTPNT 

4960  IF  XDUM(I%) *XDUM(I%-1)  THEN  4970  ELSE  4990 

4970  PRINT  "DUPLICATE  SAMPLE— PDF  ABORTS" 


4980  PDFFLAG$="S" : PDF=0 :  RETURN 

4990  PDFMIN=(G-GMINUS)/(XDUM(IS)-XDUM(IS-1)) 

5000  XLO=.5*(XDOM(I%-1)+XDOM(I%) ) 

5010  GOTO  5120 

5020  'ELSE  Segment — Interpolate  To  Next  Derivative 
5030  PDFMIN= (GPLUS-G) / (XDDM ( IS+1) -XDUM ( 1%) ) 

5040  XLO=.5* (XDDM(I%) +XDUM(I%+1) ) 

5050  IF  IS=SAMPSIZES-1  THEN  PDFMAX=0 : XHI=DUMAX :GOTO  5120 
5060  'CALL  PLOTPNT(I%+2, I PLOTS, SAMPS IZE%;PLOTPNT) 

5070  IPNTS-IS+2 

5080  GOSUB  3060 

5090  GHIGH=PLOTPNT 

5100  PDFMAX= (GHIGH-GPLUS) / (XDUM(I%+2) -XDUM(I%+1) ) 

5110  XHI=.5* (XDOM(I%+l) +XDUM(IS+2) ) 

5120  PDF=PDFMIN+ (PDFMAX-PDFMIN) / (XHI-XLO) * (X-XLO) 

5130  IF  PDF<0  THEN  PRINT  "PDF  ABORT  IN  POINTCDF" : STOP 
5140  RETURN 
5150  END 

5160  '  SUBROUTINE  DUMPDATA(XDATA[] , SIZES) 

5170  ON  ERROR  GOTO  5410 

5180  IF  MID$ (FILNAM$, 2 ,1) "  THEN  FILNAM$=MID$ (FILNAM$,3) 
5190  FILE $=  PREFIX $+  F ILNAM $ 

5200  OLDF IL  $=F ILE  $+ " . BAK * 

5210  KILL  OLDFIL$ 

5220  NAME  FILE$  AS  OLDFIL$ 

5230  OPEN  "0",1,FILE$ 

5240  NELS=5 

5250  NRECS*SIZES\NELS 

5260  IF  SIZES  MOD  NELSOO  THEN  NRECS-NRECS+l 

5270  XINDEXS=0 

5280  PRINT  #1 , SIZES ;NELS 

5290  FOR  RECORDS =1  TO  NRECS 

5300  FOR  ELEMNTS=*1  TO  NELS 

5310  X INDEXS =X INDEXS +1 

5320  IF  XINDEXS>SIZES  THEN  PRINT  #l,:CLOSE  1 : RETURN 
5330  PRINT  # 1,XDATA(X INDEXS ) ; 

5340  NEXT  ELEMNTS 
5350  PRINT  #1, 

5360  NEXT  RECORDS 
5370  CLOSE  1 
5380  RETURN 
5390  END 

5400  "SUBROUTINE  DISKERR"""""""" ' 

5410  IF  ERR=53  AND  ERL01300  THEN  RESUME  NEXT 
5420  ON  ERROR  GOTO  0 
5430  RETURN 
5440  END 

5450  '"""""SUBROUTINE  PLOTPCDF"""""""""" 

5460  MEMORY=FRE{0) 

5470  PRINT  "MEMORY  LEFT  IS  "; MEMORY; "BYTES" 

5480  MEM«1000 

5490  IF  MEMORY<MEM  THEN  5520 


5500  INPUT  "Do  You  Want  Plotdata  Written  To  Disk";ANS$ 

5510  IF  ANS$="Y"  THEN  5520  ELSE  5630 
5520  COMMON  MAXPDF,XMIN,XMAX, FILES 

5530  DIM  X(l) ,Y(1) ,Y2(1) •  Sets  Starting  Addresses 
5540  CLOSE 

5550  INPUT  "MUST  WRITE  TO  DISK.  INPUT  FILNAM$ : " , FILNAM$ 

5560  IF  MID$(FILNAM$,2, 1) =  "  THEN 

FILNAM$=MID$ (FILNAM$ r 3 ) 

5570  FILE$=PREFIX$+FILNAM$ 

5580  OLDFIL$=FILE$+" . BAK" 

5590  KILL  OLDFILS 

5600  NAME  FILE$  AS  OLDFIL$ 

5610  OPEN  "0",1,FILE$ 

5620  GOTO  5660 

5630  COMMON  X()  ,Y()  ,Y2() 

5640  FILE$="As"'  Sets  Address 

5650  DIM  X (101) , Y (101) , Y2 (101) 

5660  PRINT  "XMAX  AND  XMIN  ARE  "; XMAX; XMIN: INPUT  "CHANGE 
VALUES"; CHANGES 

5670  IF  CHANGE$="Y"  THEN  INPUT  "INPUT  XMIN, XMAX: " , XMIN, XMAX 

5680  HX= (XMAX-XMIN) /100 

5690  FOR  POINT%=l  TO  101 

5700  INDEX%=POINT%-l 

5710  X=XMIN+INDEX%*HX 

5720  GOSUB  3770 

5730  PRINT  "AT  X=";X? "CDF, PDF  ARE  " ; AVGCDF ; AVGPDF 
5740  IF  AVGPDF >MAXPDF  THEN  MAXPDF=AVGPDF 

5750  IF  MEMORY<MEM  OR  ANS$»"Y"  THEN  5760  ELSE  5780 

5760  PRINT  #1,X; AVGPDF; AVGCDF 

5770  GOTO  5790 

5780  X ( INDEX% ) =X : Y ( INDEX% ) “AVGPDF : 

Y2 ( INDEX% ) =AVGCDF 

5790  NEXT  POINT% 

5800  IF  MEMORY<MEM  OR  ANS$="Y"  THEN  5810  ELSE  5870 
5810  PRINT  #1,MAXPDF; XMIN; XMAX 

5820  CLOSE 

5830  PRINT  "MAXPDF ; XMIN ; XMAX ; F ILE$= " ; MAXPDF ; XMIN ; 

XMAX; FILES 

5840  INPUT  "Chain  In  PLOTCDF  (Y/N) ";ANS$ 

5850  IF  ANS$="N"  THEN  PRINT  "NORMAL  TERMINATION" : STOP 
5860  CHAIN  "PLOTCDF" 

5870  PRINT  "MAXPDF  IS  "; MAXPDF 

5880  INPUT  "CHANGE  MAXPDF  ";PDF$ 

5890  IF  PDF$="Y"  THEN  INPUT  "Input  MAXPDF :", MAXPDF 
5900  FOR  INDEX%=0  TO  100 

5910  Y(INDEX%)=Y(INDEX%) /MAXPDF 

5920  NEXT  INDEX% 


5930  PRINT  "MAXPDF, XMIN, XMAX  ARE 

" ; MAXPDF ; XMIN ; XMAX ; " PAUSING " ; STOP 
5940  CHAIN  "MXPLOT.002" 

5950  RETURN 

5960  END 


B.14 


5970  'SUBROUTINE  EXTRAP (SUBSAMP [] ,TOGGLE$;DELTAL,DELTAU) 
5980  ALFAPLOT=-.5:BETAPLOT=0 
5990  IF  TOGGLE$="R"  OR  PDFTOGGLE$="S"  THEN 
DELTAL=1 :  DELTAU=1 :  RETURN 

6000  IF  EXTRAPTOGGLE$="S"  THEN  DELTAL=1 :DELTAU=1 : RETURN 
6010  RLl = 2 * SUBSAMP ( 2 , SAMPLE% ) - SUBSAMP ( 1 f  SAMPLE% ) 

-SUBSAMP ( 3 , SAMPLE% ) 

6020  RL2= SUBSAMP ( 2 , SAMPLE! ) - SUBSAMP ( 1 , SAMPLE! ) 

6030  RL3* SUBSAMP ( 3 , SAMPLE! } -SUBSAMP ( 2 , SAMPLE! ) 

6040  RL4= SUB SAMP ( 3 , SAMPLE! ) - SUBSAMP ( 1 ,  SAMPLE! ) 

6050  RL=RL1*RL2/(RL3*RL4) 

6060  IF  RL>=1  THEN  6070  ELSE  6150 

6070  PRINT  "EXACT  EXTRAPOLATION  FAILS  FOR  DELTAL" : STOP 

6080  RL=RL4/RL3 

6090  IF  RL>=4  THEN  6100  ELSE  6130 

6100  PRINT  "METHOD  2  FAILS  FOR  DELTAL" sSTOP 

6110  DELTAL=2* (1+ALFAPLOT) 

6120  GOTO  6160 

6130  DELTAL=2* (1+ALFAPLOT) /( 2-. 5*RL) 

6140  GOTO  6160 

6150  DELTAL=2* ( 1+ALFAPLOT) /(1-RL) 

6160  RH1A=SUBSAMP(SAMPSIZE! , SAMPLE!) 

-SUBSAMP ( SAMPSIZ  E!-l , SAMPLE! ) 

6170  RH1B=SUBSAMP(SAMPSIZE!-1, SAMPLE!) 

- SUBSAMP ( SAMPSIZ E! -2 * SAMPLE! ) 

6180  RH1=RH1A/RH1B 

6190  RH2 A=  2 ‘SUBSAMP ( SAMPS I Z  E! -1 , SAMPLE! ) 

- SUBSAMP ( SAMPS I Z  E!  -  2 ,  SAMPLE! ) 

-SUBSAMP { SAMPSIZ  E! , SAMPLE! ) 

6200  RH2B=SUBSAMP(SAMPSIZE!, SAMPLE!) 

- SUBSAMP ( SAMPS I Z  E! - 2 , SAMPLE! ) 

6210  RH2=RH2A/RH2B 
6220  RH=RH1+RH2 

6230  IF  RH>=2  THEN  6240  ELSE  6320 

6240  PRINT  "EXACT  EXTRAPOLATION  FAILS  FOR  DELTAU":STOP 

6250  RH=RH2B/RH1B 

6260  IF  RH>=4  THEN  6270  ELSE  6300 

6270  PRINT  "METHOD  2  FAILS  FOR  DELTAU":STOP 

6280  DELTAU=2* (BETAPLOT-ALFAPLOT) 

6290  GOTO  6330 

6300  DELTAU=  2  * ( BETAPLOT-ALFAPLOT) / ( 2- . 5*RH) 

6310  GOTO  6330 

6320  DELTAU=2* (BETAPLOT-ALFAPLOT) / ( 2-RH) 

6330  RETURN 
6340  END 


Appendix  C ;  Distribution  Of  A  Pane t ion 
jonotonip  In  N  Rendon  V.ar  iafrj.es 


Overview 

This  Appendix  provides  s  d e non s t r e t  i  on  of  s  powerful 
alternative  to  random  Monte  Carlo  sampling  when  nsed  to  find 
the  distribution  of  a  function  of  one  or  more  random 
variables.  Appendix  A  provided  an  introduction  to  Sweeder's 
method  [14]  of  non-parame tr ic  estimation  of  distribution  and 
density  functions.  In  that  Appendix,  the  concept  of  a 
stylized  sample  was  presented,  and  improvements  were  made  to 
one  of  Sweeder's  basic  numerical  techniques.  These  ideas 
are  briefly  reviewed.  The  problem  of  finding  the 
distribution  of  a  function  of  a  single  random  variable  is 
considered.  It  is  shown  that  a  set  of  stylized  points  from 
the  population  of  the  random  variable  X  maps  into  a  set  of 
stylized  points  from  the  population  of  Z.  where  Z»g(X),  and 
the  function  g(x)  is  monotonic.  Two  examples  are  presented. 
Functions  of  two  independent  random  variables  are  then 
considered,  and  two  examples  again  presented.  Functions 
monotonic  in  N  random  variables  are  discussed. 

A  Review  of  Sweeder ' s  Method 

Given  a  set  of  observations  from  the  population  of  the 
random  variable  Z  with  distribution  function  F^tz),  one  can 


order  the  set  from  smallest  to  largest.  Given  this  set 
denoted  by: 


[z i) , i»l ,  m 


where  z^<z.+^  then  one  may  construct  the  sample  distribution 
function  F«(z)  as  follows: 


Fg(z)=0  ▼  z  <  z  q 


(C.l) 


Fg(  z )  “Gj+tGj+j-Gj )  (  z-z  j )  /  ( z  i  +  i-z  j )  ▼  zAlz<zi  +  1  (C.2) 


Fo(z)=l  ▼  z^z 


( C .  3  ) 


Sweeder  [14]  showed  that  such  a  sample  distribution 
function  converges  uniformly  to  the  underlying  distribution 
function  Fg ( z ) . 

In  Equation  C.2,  G  i  is  a  n o n p a r a m e t r  i  c  plotting 
position,  given  by  some  rule  such  as: 


Gi=(  i  +  a)  /  ( m+p  )  ;  -li«<Pil  (C.4) 

In  Equations  C.l  and  C.3,  Zq  and  zm  +  1  are  extrapolated 


endpoints  such  that: 

FS(z0)=G0=0 
FS(zm  +  l)=Gm  +  l=<1 

The  simplist  linear  transformation  is  given  by: 


(C.5 ) 


(C.6) 


Z=X  ( C . 7  ) 

If  X  is  a  random  variable  with  distribution  function 


Fj(x),  then  the  distribution  function  of  Z  is  given  by: 


Fjr  (z(z)  )*Fj  (  X  ) 


(C.8) 


Now,  provided  Fj(x)  is  known,  one  can  collect  the 
particular  set  of  m  values: 


{xi};  i=l , m 

in  such  a  way  that  the  x?  satisfy: 


Gi-FX<‘i>=JxofX(l)dX 


(C.9) 


where 

Fx(x0)=0  (C.10) 

By  Equations  C.2  and  C.9  one  sees  that: 

Fs(x*)-Fx(x*)  (C.ll) 

For  this  particular  set  of  points,  the  sample 
distribution  function  Fg(x)  is  exact  at  the  drawn  data 
points,  x*.  Such  a  set  of  points  (reference  Appendix  A) 
will  be  referred  to  as  a  stylized  sample.  The  definition  is 
repeated  below: 

Definition:  The  set  of  points  {x*},  i=l,m,  is  defined 

as  a  stylized  sample  from  the  population  of  X  if,  for  every 
element  of  the  set, 

FS(xi)=FX(xi}  (C.12) 

and  *i<xi+i* 

Referring  back  to  Equation  C.7,  one  sees  that  the 

distribution  function  Fg(z)  is  also  exact  at  points  z*  given 
by 

z*  =  x*  (C.13 ) 

since 

FS(zi)s*FS(xi)“Gi  (C.14) 

Hence,  the  motivation  is  provided  to  search  for  a  new 

technique  in  finding  the  distribution  of  a  function  of  a 


* 


One  nay  now  consider  the  case  where  the  random  variable 
Z  is  an  arbitrary  monotonic  function  of  the  single  random 
variable  X.  That  is. 


Z=g(X) 


(C.15) 


where  g’(z)  has  no  zeroes  throughout  the  range  of  X. 


From  the  previous  discussion,  if  one  could  somehow 
obtain  a  stylized  sample  from  the  population  of  Z,  then  the 
sample  distribution  function  Fg(z)  is  exact  at  the  points 
[z*J,  i=l,m.  With  a  minor  restriction  on  the  plotting  rule 
of  Equation  C.4,  such  a  set  of  z's  is  found  by  simply  eval¬ 
uating  the  function  g(x)  at  the  points  {x*}. 

HYPOTHESIS:  Given  Z-g(X),  with  g(x)  monotonic,  and  a 
stylized  sample  from  the  population  of  X,  (x*),  i=l,m,  then 
the  stylized  sample  from  Z,  {z*},  i=l,m;  is  found  by  evalua¬ 
ting  g  ( x )  at  the  points  [x*J,  i  =  l,m.  This  is  true  provided 


the  x^  are  drawn  by  the  rule: 

G*(i  +  a)/(m  +  p)=F_(x.)=rXi  f_(x)dx 
1  *  1  J  *0  x 


and  0-2ar=l 


(C.16) 


(C.17) 


PEOOF:  The  equation  of  total  probability  [1,6]  can 
be  used  to  find  the  density  function  f^fz)  given  the  condi¬ 
tional  density  function  f2|j(z(x)).  That  is: 

fZ(z)=/xo+1  fIU)fZ  |x<z<*>>dx  (C.l  8) 


(X>fZ 


(C.l 8) 


where  xQ  and  xn+1  satisfy: 

Fj(xq)  =0 

FX(xm+l)=1 


( C . 19 ) 


(C.20) 


V* 


Since  the  objective  is  to  draw  or  construct  a  stylized 
sample  from  the  population  of  Z,  one  might  seek  the  point  z* 
that  satisfies  the  equation: 

Gj=FZ(Zj)  (C.21) 

That  is,  z*  is  one  of  the  members  of  a  stylized  sample  from 
the  population  of  Z. 

The  distribution  function  may  be  found  by  integrating 
the  density  function.  Equation  C.18  may  be  integrated  to 
the  point  z*.  and  the  integral  over  z  taken  inside  the 
integral  over  all  x.  Recalling  the  definition  of  the  condi¬ 
tional  distribution  function  leads  to  the  result: 

vv-Ji,"*1  ti(x,Fzix<V'11  <c-“’ 

The  z*  must  satisfy  Equation  C.22  for  every  j=l,m. 

J 

Now,  P2|l^zj^  is  lcnowtt  explicitly.  Since  z  is  a  func¬ 
tion  of  x  only,  the  density  function  of  ZjX  is  the  Dirac 


delta  function  with  parameter  {  given  by: 

f=g(x)  ( C . 23 ) 

Tha  t  is: 

fZlX(z)=i(z_f >  (C.24) 

and 

FZlXUj  )=/-i6(z-f,dz  ( C.  25 ) 

By  the  properties  of  the  Dirac  delta  function  [2,11], 
F2|j(Zj)  takes  only  two  values: 

FZ|X*zj)=0  T  zj<f  (C.26) 

F2|x(z*)=!  v  z*2.f  ( C  .  27  ) 


I?.,. 


C.5 


Equations  C.26  and  C.27  can  be  put  to  use  in  Equation 
C.22  by  transforming  the  integral  in  Equation  C.22  to  an 
integral  over  f.  From  Equation  C.23, 

dfag'(x)dx  (C.28) 
x-g-1<n  (C.29) 


The  transformed  integral  is: 

One  may  now  consider  increasing  versus  decreasing  func¬ 
tions  of  X. 

Monotonical lv  Inc  re  a  sing  Funct ions .  For  monotonical ly 


increasing  functions  of  X, 

g'(x)>Ovx  (C.31) 

For  this  case  the  integrand  in  Equation  C.30  is  always 
non-negative.  By  Equation  C.26, 

Fz|x^zj)=0  T  ^zj  (C.3 2) 

and  is  unity  elsewhere.  Using  Equation  C.32  in  C.30  reduces 

that  expression  to: 

Gj=Fz(x*)=J^jxo)fI(g"1(n)df/g'  (g-1(0)  ( C .  33  ) 

provided  that 

*j<*(*a+l)  (C.34) 

Transforming  back  to  the  variable  in  x  yields: 

G  =F_U*)=fg  1(zj  )f_(x)dx  (C.35) 

J  *  J  Jx0  * 


However , 

®j"FX(xj,“Jxo  £X(x,dx 

Since  Equation  C.34  holds  and  g(x)  is  increasing 


(C.36) 
it  follows 


that 


*  (zj,<*m+l 
fX(g"1(Zj)>0 

Comparison  of  Equations  C.3S  and  C.36  implies  that 
-1/  * 

* 

or  z j=g (ij  ) 


(C.3  7 ) 
(C.3  8) 


(C.39) 

(C.40) 


The  hypothesis  is  verified  for  increasing  functions  of 


Monotonicallv  Decreasing  Funct ions.  For  monotonical ly 


decreasing  functions  of  X, 
g  '  ( x )  <  0  r  x 


( C.41 ) 


For  this  case,  the  integrand  in  Equation  C.30  is  always 
non-positive.  Since  positive  integrands  are  preferred,  one 


can  use  the  relation: 


lg'(x)  l=-g'(x) 


(C.42) 


in  Equation  C.30.  This  results  in  the  expression: 


FZ(zJ*)="/l(xr)+l)fX(8”1<^))FZ|X(zj*)d^/lg'(8"1(^))l(C'43) 

Reversing  the  limits  of  integration  leads  to: 

Exploiting  again  Equations  C.26  and  C.27: 

G =F -<*?>-(*{.  \fT(g-*(£)>d£/lg'(g_1(C)) I  ( C. 45 ) 

J  *  J  JgUatW  * 


provided  that 

Zj  <g(x0) 


(C.46) 


A  variable  transformation  back  to  x  space  yields: 

G  -F_(z*)=-fg’1<Zj  )f_(x)  Ig'(x)  |dx/lg'(x)  I  ( C . 47 ) 

j  Z  j  Jxm+l  X 


Exploiting  Equation  C.29  and  again  reversing  limits 
leads  to: 

•  /*,  . 

(C.48) 


Equation  C.48  is  just  the  statement  that: 

• .  _  _  -1.  e. 


6j-F2(zJ)-Pr  CX>g-1<zJ)}-l-Fx(g_1(z*>) 


(C.49) 

Since  Equation  C.46  holds  and  g(x)  is  decreasing,  it 

follows  that 
-1 


(Zj) >xQ 


fX(g~1(Zj))>° 


( C . 50 ) 
(C.51) 


So  there  ainst  exist  an  x j ,  such  that 


xj'>x0 


and 


or 


Zj=g(Xj , ) 


(C.52) 

(C.53) 

(C.54) 


The  question  is.  is  this  x j  .  one  of  the  members  of  the 


set  of  stylized  points  drawn  from  the  population  of  X?  It 
is  if  it  satisfies  Equation  C.36  for  j=j'.  However,  it  must 
also  satisfy  the  equation  just  derived.  Equation  C.49.  Evi¬ 
dently.  one  must  require  that 

(C.55) 

Using  Equation  C.4  in  the  above,  j  and  j'  must  satisfy 
j»m+p-2a-j'  (C.56) 

But  the  restriction  of  Equation  C.17  reduces  the  above 


Gj=»l-Gj,  for  some  j=l,m;  some  j'=l,m 


to : 


j-m  +  l-j*  (C.57) 

Examination  of  the  above  equation  (see  Table  C.l)  shows 


that  under  the  restriction  of  Equation  C.17,  the  set  of  G.'s 


reaiii  unchanged  under  the  tranaf ornation  of  Equation  C.SS. 
It  is  interesting  to  note  that  of  all  the  plotting  rules 
considered  (Table  C.2),  only  the  EDF  fails  to  satisfy  the 
condition  jj-2a*l. 

In  conclusion  then.  Equations  C.S4  and  C.57  imply  that: 
Zj“g<Xj,)“g(x»+l-j)  ( C.5  8 ) 

The  hypothesis  is  therefore  verified  for  decreasing 
functions  of  X,  completing  the  proof. 

Sample  Calculat ions.  To  illustrate  the  above  consider¬ 
ations,  two  sample  problems  will  be  worked — one  involving  a 
known  increasing  function  of  X,  the  other  a  known  decreasing 
function  of  X. 

Distribution  o  f  4  Linear  Trans  form.  An  example  of 
an  increasing  function  of  X  is: 

Z»(lnX-ax)/0x  (C.59) 

where  X  is  lognormally  distributed  with  location  parameter 
ax  and  scale  parameter  0X.  The  input  distribution  is  from 
Figure  A. 22.  The  stylized  sample  is  formed  by  finding 

Zj=  ( lnx£-<*x) /0X  v  i»l,m  (C.60) 


The  distribution  of  Z,  which  should  be  normal  with  mean  zero 
and  standard  deviation  1,  is  found  using  the  method  previ¬ 
ously  described.  The  results  are  illustrated  in  Figure  C.l. 


TABLE  C.2 


rr-, 


where  ms  usual  the  y ^  satisfy 

Gi=(i+a)/(«+p)-FI(y*)-J^^fT(y)dy  (C.67) 

and  the  restriction  of  Equation  C.17  holds.  As  usual,  y q 


and  yB+1  are  extrapolated  endpoints  such  that: 

Fy(y0)s*0  ( C .  6  8 ) 

FT(ym+l>“1  <C*69> 

Selecting  a  particular  y*  one  can  in  principle  con- 


( C .  69 ) 


struct  the  conditional  density  function  f2|y?(z) 


equation. 


fjj|yj  (z)=fj(x(i)  )  |dx(z)/dz| 


from  the 


(C.7  0) 


If  Equation  C.70  is  analytically  tractable,  the 
conditional  distribution  function  may  be  found  by  direct 


integration,  i.e. 


r2ly;(z)a,/-<DIzlyivz,az'  ^•/1' 

If  Equation  C.70  is  not  easily  integrable,  one  can 


f-1  (a')dx' 


( C .  7 1 ) 


define  the  conditional  random  variables  Z*  as: 
Zi“Zly*;  i=l,m 


(C.7  2) 


Since  there  are  m  elements  in  the  set  of  stylized  points 
from  the  population  of  T,  m  conditional  random  variables  may 
be  defined.  From  the  argument  just  completed  for  functions 
of  a  single  random  variable,  stylized  sets  of  points  for 

4  | 

each  of  the  Z  may  be  drawn  by  holding  y^  fixed  and  forming 


the  sets  from  the  rule: 

{Sj*}“(g(x? ,y*)J ;  j«l,m 


(C.73) 


variables  with  a  eleaeats  each,  requiring,  at  aost,  a-4  calls 
to  the  function  g(x,y). 

In  any  case,  the  conditional  distribution  function 
F2|y?(z)  is  readily  available. 

The  density  function  of  Z  is  given  by  Equation  C.18, 


and  the  distribution  function  by: 

Fz,‘,-fn'‘tr<y>FzlT,’>dy  <c'74> 

Since  the  function  Fy(y)  is  known,  as  are  the  functions 
FZlT(z)*  t  h e  integral  of  Equation  C.74  aay  be  solved  by 
Mellin  transfora  [9],  also  known  as  the  graphical  aethod 
[8].  The  exact  details  of  this  using  stylized  sets  are 


described  in  Appendix  D. 


Again,  the  objective  is  to  construct  a  stylized  set  of 
points  froa  the  population  of  Z.  That  is,  one  finds  the  z* 
satisfying 

Gj  =  (  j+a)  /  (m  +  p)=Fz  (  z  *)  v  j=l,m  (C.75) 

or,  one  seeks  solutions  to  the  inverse  equation, 

z*=Fz1(Gj)  (C’76) 

The  approach  now  is  direct  solution  of  Equation  C.76  by 
nuaorical  iterative  aethods.  Two  exaaples  using  this 


formulation  are  presented  at  this  tiae. 

Distribution  o_f  A  Sum  o£  IaSAI«l •  An  example  of  a 
function  monotonic  in  two  independent  random  variables  is 

Z«X2+I2  (C.77) 


where  both  X  and  T  are  independent  and  distributed  uniformly 
on  the  interval  [0,1].  Again,  the  input  distribution  for 


either  X  or  T  is  that  shovn  in  Figure  A. 18.  The  stylized 
staple  froa  the  population  of  Z  is  found  by  solution  of 
Equation  C.7 6. 

Some  tedious  analytic  work  yields  the  exact 
distribution  and  density  functions  given  respectively  by: 

Fz(z)=0  ▼  zi.0  ( C .  78) 

Fz (*)=nz/4  ▼  Oi.zi.1  (C.79) 

For  the  region  l^.z^2, 

Fz(z)=z{Sin-1[l/z1/2]-Sin_1[ ( ( z-1 ) / z ) 1/2 ] )/2 

+(z-l)1/2  (C.80) 


Fz(z)=l  ▼  z2.2  (C.81) 

fz(z)=0  v  zi.0  (C.82) 

fz(z)=n/4  v  Olzll  (C.83) 

For  the  region  l£z£2, 

fz(z)={Sin"1ll/z1/2]-Sin~1[((z-l)/z>1/2)}/2  (C.84) 

fz(z)=0  ▼  z2  2  (C.85) 


Equation  C.76  was  solved  iteratively  for  various  num¬ 
bers  of  stylized  points  by  the  method  just  described.  The 
resulting  convergence  with  increasing  number  of  stylized 
points  is  illustrated  in  Figures  C.3  through  C.6.  These 
figures  illustrate  the  overall  convergence  of  both  the  dis¬ 
tribution  and  density  functions  as  the  number  of  stylized 
points  increases.  For  purposes  of  illustration,  Sweeder's 
trigonometric  interpolation  was  used  for  these  plots.  The 
most  interesting  of  these  is  Figure  C.3  which  was  done  using 
only  S  stylized  points  from  each  of  the  input  distributions. 


Lzed  Po 


The  theory  just  developed  predicts  that  the  numerical 
a pp r o z i a a t i o a  to  the  true  d i s t r ib a t i on  function  should  be 
pinned  to  the  exact  distribution  function  at  the  percentile 
points  .10(.20)  .90.  The  crossover  points  are  very  close  to 
these,  even  though  it  is  expected  that  integration  errors 
for  such  a  saall  nuaber  of  points  Bight  lead  to  erroneous 
values  for  the  z?.  The  probability  density  function  is 
uniform  out  to  z  =  l,  and  then  has  a  longer  decaying  tail. 
The  extrapolated  endpoints  coapare  favorably  to  the  exact 
results  of  0  and  2. 

Distribution  of  A  Ratio.  Finally,  consider  the  distri¬ 
bution  of 

Z=X/Y  (C.8  6) 

where  both  X  and  Y  are  distributed  noraally  with  aean  0  and 
variance  1.  It  is  known  that  the  distribution  of  Z  is 
Cauchy  with  density  given  by: 

fz(z)=l/(n(l+z2)}  ( C . 87 ) 

and  distribution  function  given  by: 

Fz(z)=.5+Tan_1(z)/n  (C.88) 

This  was  investigated  since  the  Cauchy  is  particularly 

troublesoae  for  aoaent  propagation  aethods,  including 

Shannon  aaximua  entropy  [13],  and  because  a  hard  test  of  the 
method  was  desired  for  the  case  when  the  distribution  has 
infinite  support  and  heavy  tails.  The  results  are 
displayed  in  Figure  C.7. 


rewritten  as  combinations  of  pairs  of  independent  random 


▼arisbles.  One  then  nses  the  methods  of  the  last  section  to 
solve  a  sequence  of  binary  problems.  A  relatively  simple 
example  is: 

Z-{  (R-S)¥+T}/U  (C.89) 

where  R,S,T,U,  and  f  are  all  independently  distributed 
random  variables  with  known  distributions.  The  problem  can 
be  solved  by  solving  a  sequence  of  problems  of  two  random 
variables.  First,  one  finds  the  distribution  of 

V-R-S  (C.90) 

This  is  done  using  the  methods  just  described,  and  once  the 

distribution  of  Y  is  known.  Equation  C.89  has  been  reduced 
to 

Z  =  (VW+T)/D  ( C.91 ) 

Now  one  can  proceed  by  finding  the  distribution  of  X 

where  X  is  defined  as 

X=VW  (C.92) 

The  original  problem  has  been  reduced  in  dimension  again  to 

Z-(X+T)/0  (C.93) 

Now  X  may  be  defined  by 

Y-X+T  (C.94) 


and  the  distribution  of  Y  can  be  found 


since  the 


distributions  of  X  snd  T  are  known.  The  final  problem  to 


solve  is  thus 


Z-T/O 


(C.95) 


Bnt  this  is  easily  done,  since  the  distributions  of  T  and  0 
are  known,  and  are  independent. 

The  method  is  not  limited  to  linear  functions.  The 


multiple  random  variables,  solved  using  no  additional 
mathematical  complexity  than  that  for  two  random  variables. 
Actually,  the  requirement  of  independence  might  be  relaxed, 
depending  on  the  problem  being  considered.  If  two  of  the 
variables  are  correlated,  if  the  correlation  is  known,  and 
if  those  two  can  be  isolated  from  the  other  random 
variables,  then  the  distribution  of  their  combination  can 
be  determined.  This  is  done  in  Chapter  IV,  where  the 
distribution  of  the  safety  factor  of  a  reactor  pressure 
vessel  is  determined,  the  safety  factor  being  a  function  of 


five  random  variables,  two  of  them  correlated. 


Summary 


The  primary  use  of  algorithm  NOSWET  has  been  demon¬ 


strated  in  this  Appendix.  Sweeder's  method  of  nonparametr ic 
estimation,  as  modified  (Reference  Appendix  A),  provides  a 
new  tool  for  finding  the  distribution  of  a  function  of  a 
random  variable.  After  reviewing  Sweeder's  basic  ideas. 


-  ty 


population  of  X  saps  into  a  stylized  saaple  froa  the  popula¬ 
tion  of  Z.  The  idea  of  drawing  a  stylized  saaple  to  find 
the  distribution  of  an  output  variable  leads  to  a  nuaerical 
aethod  for  finding  the  distribution  of  a  function  of  two 
randoa  variables.  Ezaaples  were  shown,  and  the  problea  of  a 
function  of  multiple  random  variables  was  discussed.  The 
technique  provides  a  new  tool  for  the  survivability  analyst, 
since  the  stress  and  strength  distributions  are  often  func¬ 
tions  of  basic  randoa  variables  whose  input  distributions 
are  known.  The  aethod  provides  a  powerful  alternative  to 
randoa  Monte  Carlo  aethods  and  propagation  of  moaents  tech¬ 


nique  s 


I 


To  The  Reliability  Interference  Integral 

The  objective  of  this  Appendix  is  to  demonstrate  the 
utility  of  stylized  sets  of  points  (reference  Appendix  A)  in 
solving  conditional  probability  integrals,  like  that  shown 
in  Equation  C.74,  reproduced  below. 

F  z  ( * )  =J"ym+ 1  f  T  Z  |T  (D.l) 

The  integral  essentially  extracts  the  expected  value  of 

a  cumulative  distribution  function  (CDF)  with  respect  to  a 
probability  density  function  (PDF).  The  CDF  and  PDF  need 
not  be  independent,  as  Equation  D.l  above  explicitly 
indicates.  The  relationship  of  the  above  to  the  reliability 
interference  integral  may  be  seen  by  considering  T  as  the 
stress  variable  s,  and  Z|Y  as  the  strength  variable  S. 
Equation  D.l  then  becomes  the  equation  for  the  failure 
probability  pf,  given  by 

P  .=  rm+l  f  .  (s)F_(s)ds  (D.2) 

I  J  SO  S  S 

A  variable  transformation  technique,  known  as  a  Hellin 
transform  [9]  may  be  used  to  calculate  the  above  integral. 
(This  is  also  known  as  the  graphical  method  by  some  other 
authors  [8].)  The  variables  G  and  H  are  defined  by 

c<s>-/?0 


f0(S')dS' 


(D.3 ) 


H(s)=  f  (s')ds' 
Jto  * 

Then 


(D.4) 


H=f #  (s)d* 

and  Equation  D.2  may  be  written  as 


■h 


G(H)dH 


(D.5 ) 


(D.6) 


The  dependence  of  G  on  H  is  seen  through  the  relations 

Hi  =  Fa<si)  (D.7) 

or  si=F«1(Hi>  (D.8) 

At  the  particular  value  s  ^ ,  G  .  is  given  by 

Gi=FS(si)  (D.9) 

or,  writing  G  as  an  explicit  function  of  H, 


Gi=FS[Fa1(Hi)]  (D.10) 

Numerically,  one  draws  a  set  of  stylized  points  from 
the  stress  random  variable  a,  and  denotes  the  set  by 
tsT},  i=l,m 

where  Sq  and  sm  +  ^  are  the  extrapolated  endpoints  found  as 
discussed  in  Appendix  A.  Because  of  the  way  the  stylized 
points  are  drawn,  the  are  given  exactly  by 

Hi=(i+a)/(m+p) ;  -li<*lP;  i-l,m  (D.ll) 


where 


H0  =  0  (D.12) 

and  Hm+1-1  (D.13) 

The  set  of  stylized  points  from  the  strength  variable 
S,  denoted  by 

( S*} ,  i-l,m 


defines  the  strength  distribution  function  F«(s)  as  defined 


calculated  over  H  by  writing 

p  f  =Jg®+lG(H)dH  (D.14) 

The  above  integral  may  be  written  as  the  snm  of  three 
separate  integrals  on  the  intervals  [Hq,H^],  [H^,Ha],  and 
[Hm>Hm^i].  Application  of  the  trapezoidal  rule  to  each  of 
these  integrals  in  turn  yields 


pf i=<6i+Go^Hl/2 

(D.15) 

(D.16) 

Pf3  =  (Gm  +  l+Gm)(1-Hm)/2 

I 

(D.17 ) 

AH=1/ (m+0) 

(D. 18) 

The  sum  of  the  above  three  terms  is  approximately  the 


value  of  the  reliability  interference  integral 


Overview 


This  Appendix  provides  additional  information  on  the 
box-beam  wing  model  discussed  in  Chapter  VI.  In  particular, 
the  following  is  provided:  (a)  the  basic  assumptions  about 
the  box-beam  and  the  scenario  for  the  vulnerability  calcula¬ 
tions,  (b)  explicit  calculations  of  the  stress  in  each  of 
the  four  members  of  the  beam,  and  (c)  the  relative  vulnera¬ 
bility  of  each  of  the  four  components  of  the  beam  in  the 
sure-kill,  median-failure,  and  sure-safe  regions. 

Assumptions 

The  following  assumptions  are  made  regarding  the  wing 
model  and  the  scenario  for  the  calculation: 

(a)  the  wing  is  assumed  to  be  designed  for  constant 
bending  stress  along  the  length  of  the  wing  with  a 
lg  design  of  10,000  psi  (6.891C7  Pascals); 

(b)  the  wing  box  is  assumed  to  be  free  to  expand  and 
bend; 

(c)  the  effects  of  structural  discontinuities  are 
neglected; 

(d)  the  aircraft  is  directly  above  the  burst  at  burst 
time,  implying  that  only  the  lower  skin  is  heated, 
and  it  is  heated  uniformly. 


Stress  Calculat ions 

There  are  four  components  to  the  box-beam,  as  shown  in 
Figure  6.9.  These  inclnde  the  lower  skin,  the  upper  skin, 
the  lower  stringers,  and  the  upper  stringers.  Even  though 
the  lower  skin  is  the  only  heated  part,  the  effects  of  this 
heating  will  be  felt  by  all  the  members  through  the  integral 
equations  shown  in  Equations  6.25  and  6.26.  These  integrals 
have  the  explicit  values  shown  below. 

I. 


L=h  +  Ur  <z,t)dz  (E.l) 

l  j— c-e  p 

where  c  is  as  shown  in  Figure  6.9,  and  (  is  the  half¬ 
thickness  of  the  skin.  In  general  terms,  Ij  reduces  to 

I1=2tATp  (E.2) 

Similarly,  the  integral  in  Equation  6.26  is  given  by: 


>AT  (z.t)dz 

P 


(E.3) 


2  J- c— € 

This  integral  reduces  to: 

I2—2c€ATp  (E.4) 

With  these  integrals  determined,  and  denoting  the 
stress  in  the  lower  skin,  upper  skin,  lower  stringer,  and 
upper  stringer  by  a  2 »  O'  3 .  and  (T^,  respectively,  one 

finds  the  results: 

al=algML"*3582alE<V*TP 

a2— ^1  g11^  * 0 41 7  7“  1 E  ( Tp  >  ATp 

<73  =  5aigI,L/6+-5848olE(Tp)ATp 

*4—  5tflgML/6+.01519a1E(Tp)ATp 
where,  in  the  above  equations,  is  the  dimensionless  load 


(E.5) 
(E.  6 ) 
(E.7) 
<E.8) 


factor  on  the  wings,  a,  is  the  coefficient  of  linear  expan- 


sion,  E(Tp)  is  the  modulus  of  elasticity,  and  is  the  peak 
temperature  of  the  skin  at  blast  arrival  time.  Since  the 
constants  in  the  above  equations  are  dimensionless,  English 
or  metric  units  may  be  used.  For  reference  purposes,  the 
calculations  were  done  using  an  a j  of  12.7E-6  in/(in-°F)  and 
an  E(Tq)  of  1E7  psi.  (These  correspond  to  metric  values  of 
2.286E-5  m/(m-°K)  and  6.891E10  Pascals,  respectively.) 

Re lat ive  Vulnerabil itv 

The  relative  vulnerability  to  the  combined  effects  of 
both  gust  and  thermal  are  shown  in  Table  E.l.  The  three 
aircraft  are  located  in  the  sure-kill  (308  meters  range), 
m e d i a n- f a  i  1 ur e  (366.5  meters  range),  and  sure-safe  (438 
meters  range)  regions.  Each  aircraft  is  directly  above  the 
weapon  at  burst  time,  traveling  with  a  velocity  of  130 
m e t e r s  /  s e c ond  (250  knots).  The  blast  wave  catches  each 
aircraft  at  a  different  time,  resulting  in  a  gust/thermal 
synergism  that  changes  character.  The  measure  of  vulnera¬ 
bility  is  taken  as  the  point  ratio  of  the  yield  stress  of 
6061-T6  aluminum  to  the  applied  stress.  This  is  the  safety 
factor.  A  discussion  of  the  relative  vulnerabilities  in 
each  region  follows. 

The  Sure-Kil 1  Region — The  upper  skin  appears  to  be  the 
most  vulnerable  assembly  in  the  sure-kill  region.  In  each 
structure,  the  thermal  contribution  to  the  total  stress  load 


is  small,  but  this  is  particularly  true  for  the  upper  skin. 
The  upper  skin  structure  exhibits  no  pronounced  gust/thermal 
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TABLE  E.l 

RELATIVE  VULNERABILITIES  OF  A  BOX-BEAN  WING 


•>- 

REGION 

ITEM  GUST 

THERMAL 

TOTAL 

STRENGTH 

SAFETY 

LOAD 

LOAD 

STRESS 

(S-Pn) 

FACTOR 

(P») 

(Pm) 

<*-Pm) 

(S/s) 

SK 

Lover  4.13E8 
Skin 

-2 .50E7 

3 . 88E8 

2.7 1E8 

.70 

Lover  3 . 44E8 
Stringer 

4.08E7 

3 . 84E8 

2 . 7  6E8 

.72 

Upper  -3 . 44E8 

1 .06E6 

-3 . 42E8 

2 .  76E8 

.81 

•• 

Stringer 

Upper  -4.13E8 

-2 . 91E6 

-4 .  1 6E9 

2.76E8 

.66 

L_», 

Skin 

MF 

Lover  2.73E8 

-1 . 87E7 

2.56E8 

2.74E8 

1.07 

Skin 

r 

nV- 

Lover  2.30E8 

3 .05E7 

2.60E8 

2.76E8 

1.06 

V," 

Stringer 

•  *  ' 

- 

•  ►  ’ 

Upper  -2.30E8 

7 . 93E5 

-2 . 29E8 

2.7  6E8 

1.20 

N*  % 

wm*‘» 

Stringer 

u 

Upper  -2.75E8 

-2.18E6 

-2 . 77E8 

2 . 7  6E8 

1.00 

Skin 

SS 

Lover  1.99E8 
Skin 

-1.37E7 

1.8SE8 

2 . 7  6E8 

1.49 

Lover  1.67E8 
Stringer 

2 . 24E7 

1.89E8 

2 . 7  6E8 

1.46 

Upper  -1 . 67E8 

5 .82E5 

-1 .66E8 

2 . 7  6E8 

1.66 

Stringer 

Upper  -1 . 99E8 

-1 .60E6 

-2 . 01E8 

2.76E8 

1.37 

Skin 

£ 

v‘,  V 
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synergism,  since  the  load  stress  is  dominated  by  the  gust, 
and  the  strength  depends  on  fixed  room  temperature  material 
properties.  The  lower  skin  does  exhibit  gnst/thermal  syner¬ 


gism.  The  load  stress,  as  in  the  upper  skin,  is  dominated 
by  the  gust  term,  but  the  strength  properties  depend  on  the 
thermal  behavior.  Because  of  this,  this  structure  was  cho¬ 
sen  for  the  sample  calculation  done  in  Chapter  VI. 

The.  Median-Failure  RjL&Aiin. —  In  this  region,  the  most 
vulnerable  assembly  is  the  upper  skin.  This  is  because  the 
thermal  stress  contribution  acts  in  the  same  direction  (com¬ 
pressive)  as  the  gust  load  contribution.  Even  so,  the 
gust/thermal  synergism  is  slight,  the  thermal  contributing 
only  a  very  small  fraction  of  the  total  load  stress. 

In  contrast,  the  lower  stringer,  though  a  relatively 
hard  assembly,  exhibits  the  most  gust/thermal  synergism. 
The  thermal  stress  term  contributes  the  order  of  12%  of  the 
total  stress  load.  The  stress  distribution  here  would  then 
depend  on  both  the  load  factor  and  peak  temperature  random 
variables.  The  strength  function  depends  only  on  room  tem¬ 
perature  material  properties. 

The  Sur e-Sa f e  Region — Again,  the  most  vulnerable  assem¬ 
bly  is  the  upper  skin,  with  both  gust  and  thermal  stress 
terms  acting  compressively.  Again,  the  synergism  is  slight, 
since  the  thermal  contribution  remains  a  small  fraction  of 
the  total.  The  comments  made  above  regarding  the  lower 


stringer  apply  here  also 
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Chapter  III  Nomenclature 

Acronyms  And  Abbreviations 
CDFsCumulative  Distribution  Function 
DOD»Depar tment  of  Defense 

LLNLaLawrence  Livermore  National  Laboratory,  Livermore 
PDF=Probab il ity  Density  Function 
SK=A  Su’-e-Kill  Stress  Value 

yirllfrlgj 

Fg(s)=CDF  of  the  Strength  Random  Variable 
fg(s)=PDF  of  the  Strength  Random  Variable 
Fs(s)-CDF  of  the  Stress  Random  Variable 
fs(s)=PDF  of  the  Stress  Random  Variable 
P£=Failure  Probability 
S=Strength  Random  Variable 
a-Stress  Random  Variable 

s  or  S=A  particular  value  of  Stress  or  Strength 
X|«Random  Variable  Inputs 
Z^Random  Variable  Output 

i ( s-»)“Dirac  Delta  PDF  with  parameter  s 
T)=Safety  Factor  Random  Variable  (dimensionless) 
(^Margin  Random  Variable  (stress  units) 


Chanter  IV  Nomenclature 

Acronyms  And  Abbr ev iat ions 
CDF— Cumulat ive  Distribution  Function 
cascent ime  ter  s 
EPa=Eilopascals 

PDF=Probab il ity  Density  Function 
psispounds  per  square  inch 

Variables 

AaModeling  coefficient  random  variable  (dimensionless) 
a^Weibull  scale  parameter  (EPa) 
baVeibull  shape  parameter  (dimensionless) 
caWeibull  location  parameter  (EPa) 

D=Outer  vail  diameter  random  variable  (cm) 
daInner  vail  diameter  random  variable  (cm) 

Fp(p)»CDF  of  the  applied  pressure 


FCyj=Strain  functional  random  variable  (dimensionless) 


G^=ith  percentile  (Oj(.G^il) 
i=ob s erva t ion  index  integer 
m=number  of  observations 

PaApplied  pressure  random  variable  (EPa) 

P^-Burst  pressure  random  variable  (EPa) 

P0aPressure  constant  (EPa) 

p»Particular  value  of  the  applied  pressure  (EPa) 
p£aith  observation  of  the  applied  pressure  (EPa) 
t*Wall  thickness  random  variable  (cm) 

¥<*Ratio  of  diameters  random  variable  (dimensionless) 


F.2 


X^Random  Input  Variable* 

X1=Log  of  V  random  variable  (dimensionless) 

Xj^Pfodnct  of  and  ^cyl 

XjsProdnct  of  A  with  X^  (dimensionless) 

X4-Ratio  of  X3  to  P  (KPa-1) 

Z=Random  Output  Variable 

a^constant  of  the  plotting  rule  (dimensionless) 
^constant  of  the  plotting  rule  (dimensionless) 
tu=Strain  random  variable  (dimensionless) 
e^sElongation  random  variable  (dimensionless) 
T)=Safety  factor  random  variable  (dimensionless) 
0^=Ultimate  tensile  strength  random  variable  (KPa) 
(-Margin  random  variable  (KPa) 


Chanter  X  Nomenclature 


AFITsAir  Force  Institute  of  Technology*  VPAFB*  OH 
AFR^Air  Force  Regulation 

AFWL=Air  Force  Weapons  Laboratory*  Kirtland  AFB,  NM 
CDF=Cumulat ive  Distribution  Function 

DASIACaDe f ense  Atomic  Support  Information  and  Analysis  Center 
DIALOG»Tradename  of  An  Information  Retrieval  Servioe 
DNA-Defense  Nuclear  Agency 
DOD^Depar tment  of  Defense 

1 7. 

KT-Kiloton  (Energy  Equivalent  of  10  Joules) 

MC*Mission  Completion 
PDF*Probabil ity  density  function 


F.3 


psiKpounds  per  square  inch 

RK-A  sure-kill  range  value  (meters) 

RS-A  sure-safe  range  value  (meters) 

SK*>A  sure— kill  stress  value 

SS-A  sure-safe  stress  value 

S/V=Survivability/Vulnerability 

WPAFB=Wr ight-Patter son  Air  Foree  Base,  Ohio 

Variables 

2 

A=Wing  area  (meters) 

2 

Af-Fuselage  area  (meters  ) 

At=Tail  area  (meters^) 

ax=Scale  parameter  of  power  function  distribution  of  X 
B=Total  bending  moment  (nt-mt<ers) 

Bj”Fuselagc  bending  moment  (nt-meters) 

Bt»Vertical  tail  bending  moment  (nt-meters) 
bx-Shape  parameter  of  power  function  distribution 
Cdf~Fuselage  drag  coefficient  (dimensionless) 

C^tsVertical  tail  drag  coefficient  (dimensionless) 
(^•■Coefficient  of  lift  (dimensionless) 

Cj^Q^Init ial  coefficient  of  lift 

d-z  coordinate  of  aircraft  with  respect  to  burst-plane  (m) 

Z 

Fj( z ) “Cumul a t ive  Distribution  Function  of  X 
f j( z ) -Probab il i ty  Density  Function  of  X 
Lf-Fuselage  lever  arm  (meters) 

L^-Vertical  tail  lever  arm  (meters) 


n> -Load  factor  (g's — multiples  of  weight) 


P=-Randoa  Overpressure  (Pascals)  Froa  A  Nuclear  Weapon 
p=Particular  Overpressure  Value  (Pascals) 
pasAmbient  ataospherlc  pressure  (pascals) 

P£aprobability  of  failure 

r,sr“Slant  range  (aeters)  froa  a  nuclear  burst 
SsGeneral  strength  randoa  variable 
smGeneral  stress  randoa  variable 
s»Particulsr  Value  of  S  or  a 
s0~Velocity  of  sound  in  air  (aeters/sec) 
sr,r=Slant  range  (aeters)  froa  a  nuclear  burst 
S^j^sDltiaate  Load  Stress  (A  Design  Paraaeter) 
u^Peak  wind  velocity  behind  the  shock  front  (a/sec) 
nlax  y  z =Compone nt s  of  u^  in  aircraft  reference  fraae 
up=Coaponent  of  wind  perpendicular  to  vertical  tail  (a/sec 
°rx  y  velocity  components  in  aircraft  frame 

v=*Lift  wind  velocity  (aeters/ second) 

VQsInitial  aircraft  velocity  (aeters/ second) 

W=Weight  (newtons) 

XsRandoa  Ratio  of  Failure  Stress  to  Ultimate  Load  Stress 
x^Particular  Value  of  the  Randoa  Variable  X 
x-0=*x  coordinate  of  bomb  in  burst-plane  (meters) 
xt«x  coordinate  of  aircraft  in  burst-plane  (aeters) 
l»Log  safety  factor  random  variable 
y*particular  value  of  Y 

y^y  coordinate  of  boab  in  burst-plane  (meters) 
z~Standard  Noraal  Variate 


a«Angle  of  attack  (radians) 
aQsInitial  angle  of  attack  (radians) 

ax>Location  Parameter  of  Lognoraal  Distribution  of  X 
Px“Scale  parameter  of  lognoraal  distribution  of  X 
<t>(z)»CDF  of  standard  noraal  variate  evaluated  at  a 
p= Air  density  behind  the  shock  front  (kg/a^) 

“Ambient  air  density  (kg/m^) 

(i.,  =ae an  value  of  X 
I7^°standard  deviation  of  X 
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Acronyms  And  Abbreviations 


AFWL=Air  Force  weapons  Laboratory 
CDF=Cumul at ive  Distribution  Function 
KT“Kiloton  (1012  joules  energy  equivalent) 
PDF=Probability  Density  Function 
SS=Sure-Saf e 


SK“Sur e-Kil 1 


Variables 


A.t=Area  Cross  Section  of  Wing  (meter2) 
b“Intercept  of  the  SPUTTER  vs  FLUX  Regression  Line 
bv»Chordvise  Diaension  of  Ving  Section 
CTsSpecific  Heat  ( j oule s/ (kg- 0  Kelvin)) 
csHeight  of  Ving  Section  (aeters) 

E (T) -Temper a tur e  Dependent  Elastic  Modulus  (pascals) 


grsRelative  Radiated  Power  Random  Variable 


F_(a)«CDF  of  Randoa  Variable  a 


<H  <H  <H 


f ( t)*Time  Dependent  Fractional  Radiated  Power 
^"Sample  Value  of  Relative  Radiated  Power 

g(T)"Ratio  of  Yield  Stress  to  Room  Temperature  Yield  Stress 

h*  Heat  Transfer  Coefficient  (  j  onl  e  s /• -Re  lw  in-m  2  ) 

Ix"Area  Moment  of  Inertia  of  Ring  Section  (meter*) 

m=Slope  of  the  SPUTTER  vs  FLUX  Regression  Line 

N^Load  Factor  Random  Variable 

n"Number  of  Observations  of  SPUTTER 

P(r)»Time  Dependent  Radiated  Power  (watts) 

P  (W)“Yield  Dependent  Peak  Radiated  Power  (watts) 

Pj(r)-Rang«  Dependent  Failure  Probability 

E»Mean  Residual 

r"Slant  Range  (meters) 

(■Generic  Stress  Random  Variable 

T-Temperature  ("  Kelvin) 

TgK“A  Sure-Kill  Temperature  (a  Kelvin) 

TSS«A  Sure-Safe  Temperature  (°  Kelvin) 

Tp=Peak  Temperature  Random  Variable  (°-Kelvin) 

Tr“Transmi t tance  Factor 

TQ»Ambient  Temperature  (*-Kelvin) 

■■Ratio  of  Real  Temperature  to  Ambient  Temperature 

"Characteristic  Temperature 
c 

^Geometric  Part  of  Characteristic  Temperature 
t"Time  (sec) 
t  "Cooling  Time  (sec) 

w 

t  "Time  to  Peak  Radiated  Power  (sec) 


W=*Weapon  Yield  (joules) 

Yf-Value  of  FLUX  (vatti) 

Y#=*Radiated  Power  Randos  Variable  (watts) 

Ys-Mediau  Value  of  SPUTTER  (watts) 

Y^’oTheraal  Yield  (watts) 
y^Observed  SPUTTER  value  (watts) 

Z~Standard  Normal  Randoa  Variable 

z**Coordinate  (aeters)  Measured  From  Center  of  Wing  Section 
z£kA  Variate  Froa  the  Standard  Normal  Distribution 
a^Absorptivity  of  Shin  Surface 
<»j=Coe f f ic ient  of  Linear  Expansion  (°  R-1) 

OpsLocation  Parameter  of  Distribution  of  Radiated  Power 
Pp=>Scale  Parameter  of  Distribution  of  Radiated  Power 
dT=Teaperature  Difference  (#-Relvin) 

Ax=Skin  Thickness  (meters) 

C^Ving  Skin  Half  Thickness  (aeters) 

?=Ratio  of  tM  to  tc 

0(t)=Time  Dependent  Look  Angle  (Radians) 
p=>Mass  Density  (kg/m^) 

ds,F=Conditioaal  Standard  Deviation  of  SPUTTER  given  FLUX 

tf^g»Detign  Stress  (pascals) 

distress  (pascals)  Due  to  Oust  Loads 

<rgg»Sure-Eill  Stress  (pascals) 

dtll£“Stress  (pascals)  Due  to  Theraal  Loads 

distress  (pascals) 

T^Diaensionless  Time  (Ratio  of  t  to  t_ ) 
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